Changes in src/tesselation.cpp [125b3c:e359a8]
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
src/tesselation.cpp
r125b3c re359a8 9 9 10 10 #include "helpers.hpp" 11 #include "info.hpp"12 11 #include "linkedcell.hpp" 13 12 #include "log.hpp" … … 23 22 /** Constructor of BoundaryPointSet. 24 23 */ 25 BoundaryPointSet::BoundaryPointSet() : 26 LinesCount(0), 27 value(0.), 28 Nr(-1) 29 { 30 Info FunctionInfo(__func__); 31 Log() << Verbose(1) << "Adding noname." << endl; 24 BoundaryPointSet::BoundaryPointSet() 25 { 26 LinesCount = 0; 27 Nr = -1; 28 value = 0.; 32 29 }; 33 30 … … 35 32 * \param *Walker TesselPoint this boundary point represents 36 33 */ 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 Node " << *Walker << endl; 34 BoundaryPointSet::BoundaryPointSet(TesselPoint * Walker) 35 { 36 node = Walker; 37 LinesCount = 0; 38 Nr = Walker->nr; 39 value = 0.; 45 40 }; 46 41 … … 51 46 BoundaryPointSet::~BoundaryPointSet() 52 47 { 53 Info FunctionInfo(__func__); 54 //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl; 48 //Log() << Verbose(5) << "Erasing point nr. " << Nr << "." << endl; 55 49 if (!lines.empty()) 56 50 eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl; … … 63 57 void BoundaryPointSet::AddLine(class BoundaryLineSet *line) 64 58 { 65 Info FunctionInfo(__func__); 66 Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." 59 Log() << Verbose(6) << "Adding " << *this << " to line " << *line << "." 67 60 << endl; 68 61 if (line->endpoints[0] == this) … … 92 85 /** Constructor of BoundaryLineSet. 93 86 */ 94 BoundaryLineSet::BoundaryLineSet() : 95 Nr(-1) 96 { 97 Info FunctionInfo(__func__); 87 BoundaryLineSet::BoundaryLineSet() 88 { 98 89 for (int i = 0; i < 2; i++) 99 90 endpoints[i] = NULL; 91 Nr = -1; 100 92 }; 101 93 … … 107 99 BoundaryLineSet::BoundaryLineSet(class BoundaryPointSet *Point[2], const int number) 108 100 { 109 Info FunctionInfo(__func__);110 101 // set number 111 102 Nr = number; … … 115 106 Point[0]->AddLine(this); //Taken out, to check whether we can avoid unwanted double adding. 116 107 Point[1]->AddLine(this); // 117 // set skipped to false118 skipped = false;119 108 // clear triangles list 120 Log() << Verbose( 0) << "New Line with endpoints " << *this << "." << endl;109 Log() << Verbose(5) << "New Line with endpoints " << *this << "." << endl; 121 110 }; 122 111 … … 127 116 BoundaryLineSet::~BoundaryLineSet() 128 117 { 129 Info FunctionInfo(__func__);130 118 int Numbers[2]; 131 119 … … 146 134 for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++) 147 135 if ((*Runner).second == this) { 148 //Log() << Verbose( 0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;136 //Log() << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 149 137 endpoints[i]->lines.erase(Runner); 150 138 break; … … 152 140 } else { // there's just a single line left 153 141 if (endpoints[i]->lines.erase(Nr)) { 154 //Log() << Verbose( 0) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl;142 //Log() << Verbose(5) << "Removing Line Nr. " << Nr << " in boundary point " << *endpoints[i] << "." << endl; 155 143 } 156 144 } 157 145 if (endpoints[i]->lines.empty()) { 158 //Log() << Verbose( 0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl;146 //Log() << Verbose(5) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl; 159 147 if (endpoints[i] != NULL) { 160 148 delete(endpoints[i]); … … 173 161 void BoundaryLineSet::AddTriangle(class BoundaryTriangleSet *triangle) 174 162 { 175 Info FunctionInfo(__func__); 176 Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl; 163 Log() << Verbose(6) << "Add " << triangle->Nr << " to line " << *this << "." << endl; 177 164 triangles.insert(TrianglePair(triangle->Nr, triangle)); 178 165 }; … … 184 171 bool BoundaryLineSet::IsConnectedTo(class BoundaryLineSet *line) 185 172 { 186 Info FunctionInfo(__func__);187 173 if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1])) 188 174 return true; … … 199 185 bool BoundaryLineSet::CheckConvexityCriterion() 200 186 { 201 Info FunctionInfo(__func__);202 187 Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck; 203 188 // get the two triangles 204 189 if (triangles.size() != 2) { 205 eLog() << Verbose( 0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl;190 eLog() << Verbose(1) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl; 206 191 return true; 207 192 } 208 193 // check normal vectors 209 194 // have a normal vector on the base line pointing outwards 210 //Log() << Verbose( 0) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl;195 //Log() << Verbose(3) << "INFO: " << *this << " has vectors at " << *(endpoints[0]->node->node) << " and at " << *(endpoints[1]->node->node) << "." << endl; 211 196 BaseLineCenter.CopyVector(endpoints[0]->node->node); 212 197 BaseLineCenter.AddVector(endpoints[1]->node->node); … … 214 199 BaseLine.CopyVector(endpoints[0]->node->node); 215 200 BaseLine.SubtractVector(endpoints[1]->node->node); 216 //Log() << Verbose( 0) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl;201 //Log() << Verbose(3) << "INFO: Baseline is " << BaseLine << " and its center is at " << BaseLineCenter << "." << endl; 217 202 218 203 BaseLineNormal.Zero(); … … 222 207 class BoundaryPointSet *node = NULL; 223 208 for(TriangleMap::iterator runner = triangles.begin(); runner != triangles.end(); runner++) { 224 //Log() << Verbose( 0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl;209 //Log() << Verbose(3) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl; 225 210 NormalCheck.AddVector(&runner->second->NormalVector); 226 211 NormalCheck.Scale(sign); … … 229 214 BaseLineNormal.CopyVector(&runner->second->NormalVector); // yes, copy second on top of first 230 215 else { 231 eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl; 216 Log() << Verbose(1) << "CRITICAL: Triangle " << *runner->second << " has zero normal vector!" << endl; 217 exit(255); 232 218 } 233 219 node = runner->second->GetThirdEndpoint(this); 234 220 if (node != NULL) { 235 //Log() << Verbose( 0) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl;221 //Log() << Verbose(3) << "INFO: Third node for triangle " << *(runner->second) << " is " << *node << " at " << *(node->node->node) << "." << endl; 236 222 helper[i].CopyVector(node->node->node); 237 223 helper[i].SubtractVector(&BaseLineCenter); 238 224 helper[i].MakeNormalVector(&BaseLine); // we want to compare the triangle's heights' angles! 239 //Log() << Verbose( 0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl;225 //Log() << Verbose(4) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl; 240 226 i++; 241 227 } else { 242 eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl;228 //eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl; 243 229 return true; 244 230 } 245 231 } 246 //Log() << Verbose( 0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl;232 //Log() << Verbose(3) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl; 247 233 if (NormalCheck.NormSquared() < MYEPSILON) { 248 Log() << Verbose( 0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl;234 Log() << Verbose(3) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl; 249 235 return true; 250 236 } … … 252 238 double angle = GetAngle(helper[0], helper[1], BaseLineNormal); 253 239 if ((angle - M_PI) > -MYEPSILON) { 254 Log() << Verbose( 0) << "ACCEPT: Angle is greater than pi: convex." << endl;240 Log() << Verbose(3) << "ACCEPT: Angle is greater than pi: convex." << endl; 255 241 return true; 256 242 } else { 257 Log() << Verbose( 0) << "REJECT: Angle is less than pi: concave." << endl;243 Log() << Verbose(3) << "REJECT: Angle is less than pi: concave." << endl; 258 244 return false; 259 245 } … … 266 252 bool BoundaryLineSet::ContainsBoundaryPoint(class BoundaryPointSet *point) 267 253 { 268 Info FunctionInfo(__func__);269 254 for(int i=0;i<2;i++) 270 255 if (point == endpoints[i]) … … 279 264 class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(class BoundaryPointSet *point) 280 265 { 281 Info FunctionInfo(__func__);282 266 if (endpoints[0] == point) 283 267 return endpoints[1]; … … 302 286 /** Constructor for BoundaryTriangleSet. 303 287 */ 304 BoundaryTriangleSet::BoundaryTriangleSet() : 305 Nr(-1) 306 { 307 Info FunctionInfo(__func__); 288 BoundaryTriangleSet::BoundaryTriangleSet() 289 { 308 290 for (int i = 0; i < 3; i++) 309 291 { … … 311 293 lines[i] = NULL; 312 294 } 295 Nr = -1; 313 296 }; 314 297 … … 317 300 * \param number number of triangle 318 301 */ 319 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) : 320 Nr(number) 321 { 322 Info FunctionInfo(__func__); 302 BoundaryTriangleSet::BoundaryTriangleSet(class BoundaryLineSet *line[3], int number) 303 { 323 304 // set number 305 Nr = number; 324 306 // set lines 325 for (int i = 0; i < 3; i++) { 326 lines[i] = line[i]; 327 lines[i]->AddTriangle(this); 328 } 307 Log() << Verbose(5) << "New triangle " << Nr << ":" << endl; 308 for (int i = 0; i < 3; i++) 309 { 310 lines[i] = line[i]; 311 lines[i]->AddTriangle(this); 312 } 329 313 // get ascending order of endpoints 330 PointMapOrderMap;314 map<int, class BoundaryPointSet *> OrderMap; 331 315 for (int i = 0; i < 3; i++) 332 316 // for all three lines 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 } 317 for (int j = 0; j < 2; j++) 318 { // for both endpoints 319 OrderMap.insert(pair<int, class BoundaryPointSet *> ( 320 line[i]->endpoints[j]->Nr, line[i]->endpoints[j])); 321 // and we don't care whether insertion fails 322 } 338 323 // set endpoints 339 324 int Counter = 0; 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 } 325 Log() << Verbose(6) << " with end points "; 326 for (map<int, class BoundaryPointSet *>::iterator runner = OrderMap.begin(); runner 327 != OrderMap.end(); runner++) 328 { 329 endpoints[Counter] = runner->second; 330 Log() << Verbose(0) << " " << *endpoints[Counter]; 331 Counter++; 332 } 333 if (Counter < 3) 334 { 335 eLog() << Verbose(0) << "ERROR! We have a triangle with only two distinct endpoints!" << endl; 336 performCriticalExit(); 337 } 338 Log() << Verbose(0) << "." << endl; 350 339 }; 351 340 … … 356 345 BoundaryTriangleSet::~BoundaryTriangleSet() 357 346 { 358 Info FunctionInfo(__func__);359 347 for (int i = 0; i < 3; i++) { 360 348 if (lines[i] != NULL) { 361 349 if (lines[i]->triangles.erase(Nr)) { 362 //Log() << Verbose( 0) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl;350 //Log() << Verbose(5) << "Triangle Nr." << Nr << " erased in line " << *lines[i] << "." << endl; 363 351 } 364 352 if (lines[i]->triangles.empty()) { 365 //Log() << Verbose( 0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;353 //Log() << Verbose(5) << *lines[i] << " is no more attached to any triangle, erasing." << endl; 366 354 delete (lines[i]); 367 355 lines[i] = NULL; … … 369 357 } 370 358 } 371 //Log() << Verbose( 0) << "Erasing triangle Nr." << Nr << " itself." << endl;359 //Log() << Verbose(5) << "Erasing triangle Nr." << Nr << " itself." << endl; 372 360 }; 373 361 … … 378 366 void BoundaryTriangleSet::GetNormalVector(Vector &OtherVector) 379 367 { 380 Info FunctionInfo(__func__);381 368 // get normal vector 382 369 NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node); … … 385 372 if (NormalVector.ScalarProduct(&OtherVector) > 0.) 386 373 NormalVector.Scale(-1.); 387 Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl;388 374 }; 389 375 … … 402 388 bool BoundaryTriangleSet::GetIntersectionInsideTriangle(Vector *MolCenter, Vector *x, Vector *Intersection) 403 389 { 404 Info FunctionInfo(__func__);405 390 Vector CrossPoint; 406 391 Vector helper; 407 392 408 393 if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) { 409 eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl;394 Log() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl; 410 395 return false; 411 396 } … … 423 408 } while (CrossPoint.NormSquared() < MYEPSILON); 424 409 if (i==3) { 425 eLog() << Verbose(0) << "Could not find any cross points, something's utterly wrong here!" << endl; 410 eLog() << Verbose(1) << "Could not find any cross points, something's utterly wrong here!" << endl; 411 exit(255); 426 412 } 427 413 CrossPoint.SubtractVector(endpoints[i%3]->node->node); // cross point was returned as absolute vector … … 442 428 bool BoundaryTriangleSet::ContainsBoundaryLine(class BoundaryLineSet *line) 443 429 { 444 Info FunctionInfo(__func__);445 430 for(int i=0;i<3;i++) 446 431 if (line == lines[i]) … … 455 440 bool BoundaryTriangleSet::ContainsBoundaryPoint(class BoundaryPointSet *point) 456 441 { 457 Info FunctionInfo(__func__);458 442 for(int i=0;i<3;i++) 459 443 if (point == endpoints[i]) … … 468 452 bool BoundaryTriangleSet::ContainsBoundaryPoint(class TesselPoint *point) 469 453 { 470 Info FunctionInfo(__func__);471 454 for(int i=0;i<3;i++) 472 455 if (point == endpoints[i]->node) … … 481 464 bool BoundaryTriangleSet::IsPresentTupel(class BoundaryPointSet *Points[3]) 482 465 { 483 Info FunctionInfo(__func__);484 466 return (((endpoints[0] == Points[0]) 485 467 || (endpoints[0] == Points[1]) … … 503 485 bool BoundaryTriangleSet::IsPresentTupel(class BoundaryTriangleSet *T) 504 486 { 505 Info FunctionInfo(__func__);506 487 return (((endpoints[0] == T->endpoints[0]) 507 488 || (endpoints[0] == T->endpoints[1]) … … 525 506 class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(class BoundaryLineSet *line) 526 507 { 527 Info FunctionInfo(__func__);528 508 // sanity check 529 509 if (!ContainsBoundaryLine(line)) … … 542 522 void BoundaryTriangleSet::GetCenter(Vector *center) 543 523 { 544 Info FunctionInfo(__func__);545 524 center->Zero(); 546 525 for(int i=0;i<3;i++) … … 555 534 ostream &operator <<(ostream &ost, const BoundaryTriangleSet &a) 556 535 { 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 << "]"; 536 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," 537 << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]"; 560 538 return ost; 561 539 }; 562 540 563 // ======================================== Polygons on Boundary =================================564 565 /** Constructor for BoundaryPolygonSet.566 */567 BoundaryPolygonSet::BoundaryPolygonSet() :568 Nr(-1)569 {570 Info FunctionInfo(__func__);571 };572 573 /** Destructor of BoundaryPolygonSet.574 * Just clears endpoints.575 * \note When removing triangles from a class Tesselation, use RemoveTesselationTriangle()576 */577 BoundaryPolygonSet::~BoundaryPolygonSet()578 {579 Info FunctionInfo(__func__);580 endpoints.clear();581 Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl;582 };583 584 /** Calculates the normal vector for this triangle.585 * Is made unique by comparison with \a OtherVector to point in the other direction.586 * \param &OtherVector direction vector to make normal vector unique.587 * \return allocated vector in normal direction588 */589 Vector * BoundaryPolygonSet::GetNormalVector(const Vector &OtherVector) const590 {591 Info FunctionInfo(__func__);592 // get normal vector593 Vector TemporaryNormal;594 Vector *TotalNormal = new Vector;595 PointSet::const_iterator Runner[3];596 for (int i=0;i<3; i++) {597 Runner[i] = endpoints.begin();598 for (int j = 0; j<i; j++) { // go as much further599 Runner[i]++;600 if (Runner[i] == endpoints.end()) {601 eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl;602 performCriticalExit();603 }604 }605 }606 TotalNormal->Zero();607 int counter=0;608 for (; Runner[2] != endpoints.end(); ) {609 TemporaryNormal.MakeNormalVector((*Runner[0])->node->node, (*Runner[1])->node->node, (*Runner[2])->node->node);610 for (int i=0;i<3;i++) // increase each of them611 Runner[i]++;612 TotalNormal->AddVector(&TemporaryNormal);613 }614 TotalNormal->Scale(1./(double)counter);615 616 // make it always point inward (any offset vector onto plane projected onto normal vector suffices)617 if (TotalNormal->ScalarProduct(&OtherVector) > 0.)618 TotalNormal->Scale(-1.);619 Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl;620 621 return TotalNormal;622 };623 624 /** Calculates the center point of the triangle.625 * Is third of the sum of all endpoints.626 * \param *center central point on return.627 */628 void BoundaryPolygonSet::GetCenter(Vector * const center) const629 {630 Info FunctionInfo(__func__);631 center->Zero();632 int counter = 0;633 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {634 center->AddVector((*Runner)->node->node);635 counter++;636 }637 center->Scale(1./(double)counter);638 Log() << Verbose(1) << "Center is at " << *center << "." << endl;639 }640 641 /** Checks whether the polygons contains all three endpoints of the triangle.642 * \param *triangle triangle to test643 * \return true - triangle is contained polygon, false - is not644 */645 bool BoundaryPolygonSet::ContainsBoundaryTriangle(const BoundaryTriangleSet * const triangle) const646 {647 Info FunctionInfo(__func__);648 return ContainsPresentTupel(triangle->endpoints, 3);649 };650 651 /** Checks whether the polygons contains both endpoints of the line.652 * \param *line line to test653 * \return true - line is of the triangle, false - is not654 */655 bool BoundaryPolygonSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const656 {657 Info FunctionInfo(__func__);658 return ContainsPresentTupel(line->endpoints, 2);659 };660 661 /** Checks whether point is any of the three endpoints this triangle contains.662 * \param *point point to test663 * \return true - point is of the triangle, false - is not664 */665 bool BoundaryPolygonSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const666 {667 Info FunctionInfo(__func__);668 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {669 Log() << Verbose(0) << "Checking against " << **Runner << endl;670 if (point == (*Runner)) {671 Log() << Verbose(0) << " Contained." << endl;672 return true;673 }674 }675 Log() << Verbose(0) << " Not contained." << endl;676 return false;677 };678 679 /** Checks whether point is any of the three endpoints this triangle contains.680 * \param *point TesselPoint to test681 * \return true - point is of the triangle, false - is not682 */683 bool BoundaryPolygonSet::ContainsBoundaryPoint(const TesselPoint * const point) const684 {685 Info FunctionInfo(__func__);686 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)687 if (point == (*Runner)->node) {688 Log() << Verbose(0) << " Contained." << endl;689 return true;690 }691 Log() << Verbose(0) << " Not contained." << endl;692 return false;693 };694 695 /** Checks whether given array of \a *Points coincide with polygons's endpoints.696 * \param **Points pointer to an array of BoundaryPointSet697 * \param dim dimension of array698 * \return true - set of points is contained in polygon, false - is not699 */700 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const701 {702 Info FunctionInfo(__func__);703 int counter = 0;704 Log() << Verbose(1) << "Polygon is " << *this << endl;705 for(int i=0;i<dim;i++) {706 Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl;707 if (ContainsBoundaryPoint(Points[i])) {708 counter++;709 }710 }711 712 if (counter == dim)713 return true;714 else715 return false;716 };717 718 /** Checks whether given PointList coincide with polygons's endpoints.719 * \param &endpoints PointList720 * \return true - set of points is contained in polygon, false - is not721 */722 bool BoundaryPolygonSet::ContainsPresentTupel(const PointSet &endpoints) const723 {724 Info FunctionInfo(__func__);725 size_t counter = 0;726 Log() << Verbose(1) << "Polygon is " << *this << endl;727 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {728 Log() << Verbose(1) << " Testing endpoint " << **Runner << endl;729 if (ContainsBoundaryPoint(*Runner))730 counter++;731 }732 733 if (counter == endpoints.size())734 return true;735 else736 return false;737 };738 739 /** Checks whether given set of \a *Points coincide with polygons's endpoints.740 * \param *P pointer to BoundaryPolygonSet741 * \return true - is the very triangle, false - is not742 */743 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const744 {745 return ContainsPresentTupel((const PointSet)P->endpoints);746 };747 748 /** Gathers all the endpoints' triangles in a unique set.749 * \return set of all triangles750 */751 TriangleSet * BoundaryPolygonSet::GetAllContainedTrianglesFromEndpoints() const752 {753 Info FunctionInfo(__func__);754 pair <TriangleSet::iterator, bool> Tester;755 TriangleSet *triangles = new TriangleSet;756 757 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)758 for(LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)759 for(TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) {760 //Log() << Verbose(0) << " Testing triangle " << *(Sprinter->second) << endl;761 if (ContainsBoundaryTriangle(Sprinter->second)) {762 Tester = triangles->insert(Sprinter->second);763 if (Tester.second)764 Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl;765 }766 }767 768 Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl;769 return triangles;770 };771 772 /** Fills the endpoints of this polygon from the triangles attached to \a *line.773 * \param *line lines with triangles attached774 * \return true - polygon contains endpoints, false - line was NULL775 */776 bool BoundaryPolygonSet::FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line)777 {778 Info FunctionInfo(__func__);779 pair <PointSet::iterator, bool> Tester;780 if (line == NULL)781 return false;782 Log() << Verbose(1) << "Filling polygon from line " << *line << endl;783 for(TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {784 for (int i=0;i<3;i++) {785 Tester = endpoints.insert((Runner->second)->endpoints[i]);786 if (Tester.second)787 Log() << Verbose(1) << " Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl;788 }789 }790 791 return true;792 };793 794 /** output operator for BoundaryPolygonSet.795 * \param &ost output stream796 * \param &a boundary polygon797 */798 ostream &operator <<(ostream &ost, const BoundaryPolygonSet &a)799 {800 ost << "[" << a.Nr << "|";801 for(PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) {802 ost << (*Runner)->node->Name;803 Runner++;804 if (Runner != a.endpoints.end())805 ost << ",";806 }807 ost<< "]";808 return ost;809 };810 811 541 // =========================================================== class TESSELPOINT =========================================== 812 542 … … 815 545 TesselPoint::TesselPoint() 816 546 { 817 Info FunctionInfo(__func__);818 547 node = NULL; 819 548 nr = -1; … … 825 554 TesselPoint::~TesselPoint() 826 555 { 827 Info FunctionInfo(__func__);828 556 }; 829 557 … … 840 568 ostream & TesselPoint::operator << (ostream &ost) 841 569 { 842 Info FunctionInfo(__func__); 843 ost << "[" << (nr) << "|" << this << "]"; 570 ost << "[" << (Name) << "|" << this << "]"; 844 571 return ost; 845 572 }; … … 852 579 PointCloud::PointCloud() 853 580 { 854 Info FunctionInfo(__func__); 581 855 582 }; 856 583 … … 859 586 PointCloud::~PointCloud() 860 587 { 861 Info FunctionInfo(__func__); 588 862 589 }; 863 590 … … 866 593 /** Constructor of class CandidateForTesselation. 867 594 */ 868 CandidateForTesselation::CandidateForTesselation (BoundaryLineSet* line) : 869 BaseLine(line), 870 ShortestAngle(2.*M_PI), 871 OtherShortestAngle(2.*M_PI) 872 { 873 Info FunctionInfo(__func__); 874 }; 875 876 877 /** Constructor of class CandidateForTesselation. 878 */ 879 CandidateForTesselation::CandidateForTesselation (TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) : 880 BaseLine(line), 881 ShortestAngle(2.*M_PI), 882 OtherShortestAngle(2.*M_PI) 883 { 884 Info FunctionInfo(__func__); 595 CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) { 596 point = candidate; 597 BaseLine = line; 885 598 OptCenter.CopyVector(&OptCandidateCenter); 886 599 OtherOptCenter.CopyVector(&OtherOptCandidateCenter); … … 890 603 */ 891 604 CandidateForTesselation::~CandidateForTesselation() { 605 point = NULL; 892 606 BaseLine = NULL; 893 607 }; 894 608 895 /** output operator for CandidateForTesselation.896 * \param &ost output stream897 * \param &a boundary line898 */899 ostream & operator <<(ostream &ost, const CandidateForTesselation &a)900 {901 ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->Name << "," << a.BaseLine->endpoints[1]->node->Name << "] with ";902 if (a.pointlist.empty())903 ost << "no candidate.";904 else {905 ost << "candidate";906 if (a.pointlist.size() != 1)907 ost << "s ";908 else909 ost << " ";910 for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++)911 ost << *(*Runner) << " ";912 ost << " at angle " << (a.ShortestAngle)<< ".";913 }914 915 return ost;916 };917 918 919 609 // =========================================================== class TESSELATION =========================================== 920 610 921 611 /** Constructor of class Tesselation. 922 612 */ 923 Tesselation::Tesselation() : 924 PointsOnBoundaryCount(0), 925 LinesOnBoundaryCount(0), 926 TrianglesOnBoundaryCount(0), 927 LastTriangle(NULL), 928 TriangleFilesWritten(0), 929 InternalPointer(PointsOnBoundary.begin()) 930 { 931 Info FunctionInfo(__func__); 613 Tesselation::Tesselation() 614 { 615 PointsOnBoundaryCount = 0; 616 LinesOnBoundaryCount = 0; 617 TrianglesOnBoundaryCount = 0; 618 InternalPointer = PointsOnBoundary.begin(); 619 LastTriangle = NULL; 620 TriangleFilesWritten = 0; 932 621 } 933 622 ; … … 938 627 Tesselation::~Tesselation() 939 628 { 940 Info FunctionInfo(__func__); 941 Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl; 629 Log() << Verbose(1) << "Free'ing TesselStruct ... " << endl; 942 630 for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { 943 631 if (runner->second != NULL) { … … 947 635 eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl; 948 636 } 949 Log() << Verbose( 0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl;637 Log() << Verbose(1) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl; 950 638 } 951 639 ; … … 956 644 Vector * Tesselation::GetCenter(ofstream *out) const 957 645 { 958 Info FunctionInfo(__func__);959 646 Vector *Center = new Vector(0.,0.,0.); 960 647 int num=0; … … 972 659 TesselPoint * Tesselation::GetPoint() const 973 660 { 974 Info FunctionInfo(__func__);975 661 return (InternalPointer->second->node); 976 662 }; … … 981 667 TesselPoint * Tesselation::GetTerminalPoint() const 982 668 { 983 Info FunctionInfo(__func__);984 669 PointMap::const_iterator Runner = PointsOnBoundary.end(); 985 670 Runner--; … … 992 677 void Tesselation::GoToNext() const 993 678 { 994 Info FunctionInfo(__func__);995 679 if (InternalPointer != PointsOnBoundary.end()) 996 680 InternalPointer++; … … 1002 686 void Tesselation::GoToPrevious() const 1003 687 { 1004 Info FunctionInfo(__func__);1005 688 if (InternalPointer != PointsOnBoundary.begin()) 1006 689 InternalPointer--; … … 1012 695 void Tesselation::GoToFirst() const 1013 696 { 1014 Info FunctionInfo(__func__);1015 697 InternalPointer = PointsOnBoundary.begin(); 1016 698 }; … … 1021 703 void Tesselation::GoToLast() const 1022 704 { 1023 Info FunctionInfo(__func__);1024 705 InternalPointer = PointsOnBoundary.end(); 1025 706 InternalPointer--; … … 1031 712 bool Tesselation::IsEmpty() const 1032 713 { 1033 Info FunctionInfo(__func__);1034 714 return (PointsOnBoundary.empty()); 1035 715 }; … … 1040 720 bool Tesselation::IsEnd() const 1041 721 { 1042 Info FunctionInfo(__func__);1043 722 return (InternalPointer == PointsOnBoundary.end()); 1044 723 }; … … 1053 732 Tesselation::GuessStartingTriangle() 1054 733 { 1055 Info FunctionInfo(__func__);1056 734 // 4b. create a starting triangle 1057 735 // 4b1. create all distances … … 1100 778 baseline->second.first->second->node->node, 1101 779 baseline->second.second->second->node->node); 1102 Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl; 780 Log() << Verbose(2) << "Plane vector of candidate triangle is "; 781 PlaneVector.Output(); 782 Log() << Verbose(0) << endl; 1103 783 // 4. loop over all points 1104 784 double sign = 0.; … … 1116 796 if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok 1117 797 continue; 1118 Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl; 798 Log() << Verbose(3) << "Projection of " << checker->second->node->Name 799 << " yields distance of " << distance << "." << endl; 1119 800 tmp = distance / fabs(distance); 1120 801 // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle) … … 1169 850 if (checker == PointsOnBoundary.end()) 1170 851 { 1171 Log() << Verbose( 2) << "Looks like we have a candidate!" << endl;852 Log() << Verbose(0) << "Looks like we have a candidate!" << endl; 1172 853 break; 1173 854 } … … 1199 880 else 1200 881 { 1201 eLog() << Verbose(0) << "No starting triangle found." << endl; 882 Log() << Verbose(1) << "No starting triangle found." << endl; 883 exit(255); 1202 884 } 1203 885 } … … 1219 901 void Tesselation::TesselateOnBoundary(const PointCloud * const cloud) 1220 902 { 1221 Info FunctionInfo(__func__);1222 903 bool flag; 1223 904 PointMap::iterator winner; … … 1238 919 // get peak point with respect to this base line's only triangle 1239 920 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far 1240 Log() << Verbose( 0) << "Current baseline is between " << *(baseline->second) << "." << endl;921 Log() << Verbose(2) << "Current baseline is between " << *(baseline->second) << "." << endl; 1241 922 for (int i = 0; i < 3; i++) 1242 923 if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1])) 1243 924 peak = BTS->endpoints[i]; 1244 Log() << Verbose( 1) << " and has peak " << *peak << "." << endl;925 Log() << Verbose(3) << " and has peak " << *peak << "." << endl; 1245 926 1246 927 // prepare some auxiliary vectors … … 1257 938 CenterVector.AddVector(BTS->endpoints[i]->node->node); 1258 939 CenterVector.Scale(1. / 3.); 1259 Log() << Verbose( 2) << "CenterVector of base triangle is " << CenterVector << endl;940 Log() << Verbose(4) << "CenterVector of base triangle is " << CenterVector << endl; 1260 941 1261 942 // normal vector of triangle … … 1264 945 BTS->GetNormalVector(NormalVector); 1265 946 NormalVector.CopyVector(&BTS->NormalVector); 1266 Log() << Verbose( 2) << "NormalVector of base triangle is " << NormalVector << endl;947 Log() << Verbose(4) << "NormalVector of base triangle is " << NormalVector << endl; 1267 948 1268 949 // vector in propagation direction (out of triangle) … … 1271 952 TempVector.CopyVector(&CenterVector); 1272 953 TempVector.SubtractVector(baseline->second->endpoints[0]->node->node); // TempVector is vector on triangle plane pointing from one baseline egde towards center! 1273 //Log() << Verbose( 0) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl;954 //Log() << Verbose(2) << "Projection of propagation onto temp: " << PropagationVector.Projection(&TempVector) << "." << endl; 1274 955 if (PropagationVector.ScalarProduct(&TempVector) > 0) // make sure normal propagation vector points outward from baseline 1275 956 PropagationVector.Scale(-1.); 1276 Log() << Verbose( 2) << "PropagationVector of base triangle is " << PropagationVector << endl;957 Log() << Verbose(4) << "PropagationVector of base triangle is " << PropagationVector << endl; 1277 958 winner = PointsOnBoundary.end(); 1278 959 … … 1280 961 for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) { 1281 962 if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints 1282 Log() << Verbose( 1) << "Target point is " << *(target->second) << ":" << endl;963 Log() << Verbose(3) << "Target point is " << *(target->second) << ":" << endl; 1283 964 1284 965 // first check direction, so that triangles don't intersect … … 1287 968 VirtualNormalVector.ProjectOntoPlane(&NormalVector); 1288 969 TempAngle = VirtualNormalVector.Angle(&PropagationVector); 1289 Log() << Verbose( 2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl;970 Log() << Verbose(4) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl; 1290 971 if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees) 1291 Log() << Verbose( 2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl;972 Log() << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl; 1292 973 continue; 1293 974 } else 1294 Log() << Verbose( 2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl;975 Log() << Verbose(4) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl; 1295 976 1296 977 // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle) … … 1298 979 LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first); 1299 980 if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) { 1300 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;981 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; 1301 982 continue; 1302 983 } 1303 984 if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) { 1304 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;985 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; 1305 986 continue; 1306 987 } … … 1319 1000 helper.ProjectOntoPlane(&TempVector); 1320 1001 if (fabs(helper.NormSquared()) < MYEPSILON) { 1321 Log() << Verbose( 2) << "Chosen set of vectors is linear dependent." << endl;1002 Log() << Verbose(4) << "Chosen set of vectors is linear dependent." << endl; 1322 1003 continue; 1323 1004 } … … 1336 1017 // calculate angle 1337 1018 TempAngle = NormalVector.Angle(&VirtualNormalVector); 1338 Log() << Verbose( 2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl;1019 Log() << Verbose(4) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl; 1339 1020 if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner 1340 1021 SmallestAngle = TempAngle; 1341 1022 winner = target; 1342 Log() << Verbose( 2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;1023 Log() << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl; 1343 1024 } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle) 1344 1025 // hence, check the angles to some normal direction from our base line but in this common plane of both targets... … … 1358 1039 SmallestAngle = TempAngle; 1359 1040 winner = target; 1360 Log() << Verbose( 2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl;1041 Log() << Verbose(4) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl; 1361 1042 } else 1362 Log() << Verbose( 2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl;1043 Log() << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl; 1363 1044 } else 1364 Log() << Verbose( 2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl;1045 Log() << Verbose(4) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl; 1365 1046 } 1366 1047 } // end of loop over all boundary points … … 1368 1049 // 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 1369 1050 if (winner != PointsOnBoundary.end()) { 1370 Log() << Verbose( 0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl;1051 Log() << Verbose(2) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl; 1371 1052 // create the lins of not yet present 1372 1053 BLS[0] = baseline->second; … … 1398 1079 TrianglesOnBoundaryCount++; 1399 1080 } else { 1400 eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl;1081 Log() << Verbose(1) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl; 1401 1082 } 1402 1083 1403 1084 // 5d. If the set of lines is not yet empty, go to 5. and continue 1404 1085 } else 1405 Log() << Verbose( 0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl;1086 Log() << Verbose(2) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl; 1406 1087 } while (flag); 1407 1088 … … 1418 1099 bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC) 1419 1100 { 1420 Info FunctionInfo(__func__);1421 1101 Vector Intersection, Normal; 1422 1102 TesselPoint *Walker = NULL; … … 1425 1105 bool AddFlag = false; 1426 1106 LinkedCell *BoundaryPoints = NULL; 1107 1108 Log() << Verbose(1) << "Begin of InsertStraddlingPoints" << endl; 1427 1109 1428 1110 cloud->GoToFirst(); … … 1435 1117 } 1436 1118 Walker = cloud->GetPoint(); 1437 Log() << Verbose( 0) << "Current point is " << *Walker << "." << endl;1119 Log() << Verbose(2) << "Current point is " << *Walker << "." << endl; 1438 1120 // get the next triangle 1439 1121 triangles = FindClosestTrianglesToPoint(Walker->node, BoundaryPoints); 1440 1122 BTS = triangles->front(); 1441 1123 if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) { 1442 Log() << Verbose( 0) << "No triangles found, probably a tesselation point itself." << endl;1124 Log() << Verbose(2) << "No triangles found, probably a tesselation point itself." << endl; 1443 1125 cloud->GoToNext(); 1444 1126 continue; 1445 1127 } else { 1446 1128 } 1447 Log() << Verbose( 0) << "Closest triangle is " << *BTS << "." << endl;1129 Log() << Verbose(2) << "Closest triangle is " << *BTS << "." << endl; 1448 1130 // get the intersection point 1449 1131 if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) { 1450 Log() << Verbose( 0) << "We have an intersection at " << Intersection << "." << endl;1132 Log() << Verbose(2) << "We have an intersection at " << Intersection << "." << endl; 1451 1133 // we have the intersection, check whether in- or outside of boundary 1452 1134 if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) { 1453 1135 // inside, next! 1454 Log() << Verbose( 0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl;1136 Log() << Verbose(2) << *Walker << " is inside wrt triangle " << *BTS << "." << endl; 1455 1137 } else { 1456 1138 // outside! 1457 Log() << Verbose( 0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl;1139 Log() << Verbose(2) << *Walker << " is outside wrt triangle " << *BTS << "." << endl; 1458 1140 class BoundaryLineSet *OldLines[3], *NewLines[3]; 1459 1141 class BoundaryPointSet *OldPoints[3], *NewPoint; … … 1465 1147 Normal.CopyVector(&BTS->NormalVector); 1466 1148 // add Walker to boundary points 1467 Log() << Verbose( 0) << "Adding " << *Walker << " to BoundaryPoints." << endl;1149 Log() << Verbose(2) << "Adding " << *Walker << " to BoundaryPoints." << endl; 1468 1150 AddFlag = true; 1469 1151 if (AddBoundaryPoint(Walker,0)) … … 1472 1154 continue; 1473 1155 // remove triangle 1474 Log() << Verbose( 0) << "Erasing triangle " << *BTS << "." << endl;1156 Log() << Verbose(2) << "Erasing triangle " << *BTS << "." << endl; 1475 1157 TrianglesOnBoundary.erase(BTS->Nr); 1476 1158 delete(BTS); … … 1480 1162 BPS[1] = OldPoints[i]; 1481 1163 NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1482 Log() << Verbose( 1) << "Creating new line " << *NewLines[i] << "." << endl;1164 Log() << Verbose(3) << "Creating new line " << *NewLines[i] << "." << endl; 1483 1165 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one 1484 1166 LinesOnBoundaryCount++; … … 1491 1173 if (NewLines[j]->IsConnectedTo(BLS[0])) { 1492 1174 if (n>2) { 1493 eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl;1175 Log() << Verbose(1) << BLS[0] << " connects to all of the new lines?!" << endl; 1494 1176 return false; 1495 1177 } else … … 1502 1184 BTS->GetNormalVector(Normal); 1503 1185 Normal.Scale(-1.); 1504 Log() << Verbose( 0) << "Created new triangle " << *BTS << "." << endl;1186 Log() << Verbose(2) << "Created new triangle " << *BTS << "." << endl; 1505 1187 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1506 1188 TrianglesOnBoundaryCount++; … … 1516 1198 // exit 1517 1199 delete(Center); 1200 Log() << Verbose(1) << "End of InsertStraddlingPoints" << endl; 1518 1201 return true; 1519 1202 }; … … 1526 1209 bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n) 1527 1210 { 1528 Info FunctionInfo(__func__);1529 1211 PointTestPair InsertUnique; 1530 1212 BPS[n] = new class BoundaryPointSet(Walker); … … 1548 1230 void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n) 1549 1231 { 1550 Info FunctionInfo(__func__);1551 1232 PointTestPair InsertUnique; 1552 1233 TPS[n] = new class BoundaryPointSet(Candidate); … … 1556 1237 } else { 1557 1238 delete TPS[n]; 1558 Log() << Verbose( 0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl;1239 Log() << Verbose(4) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl; 1559 1240 TPS[n] = (InsertUnique.first)->second; 1560 1241 } … … 1569 1250 void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const 1570 1251 { 1571 Info FunctionInfo(__func__);1572 1252 PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr); 1573 1253 if (FindPoint != PointsOnBoundary.end()) … … 1587 1267 bool insertNewLine = true; 1588 1268 1589 LineMap::iterator FindLine = a->lines.find(b->node->nr); 1590 if (FindLine != a->lines.end()) { 1591 Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl; 1592 1269 if (a->lines.find(b->node->nr) != a->lines.end()) { 1270 LineMap::iterator FindLine = a->lines.find(b->node->nr); 1593 1271 pair<LineMap::iterator,LineMap::iterator> FindPair; 1594 1272 FindPair = a->lines.equal_range(b->node->nr); 1273 Log() << Verbose(5) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl; 1595 1274 1596 1275 for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) { … … 1598 1277 if (FindLine->second->triangles.size() < 2) { 1599 1278 insertNewLine = false; 1600 Log() << Verbose( 0) << "Using existing line " << *FindLine->second << endl;1279 Log() << Verbose(4) << "Using existing line " << *FindLine->second << endl; 1601 1280 1602 1281 BPS[0] = FindLine->second->endpoints[0]; 1603 1282 BPS[1] = FindLine->second->endpoints[1]; 1604 1283 BLS[n] = FindLine->second; 1605 1606 // remove existing line from OpenLines1607 CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]);1608 if (CandidateLine != OpenLines.end()) {1609 Log() << Verbose(1) << " Removing line from OpenLines." << endl;1610 delete(CandidateLine->second);1611 OpenLines.erase(CandidateLine);1612 } else {1613 eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl;1614 }1615 1284 1616 1285 break; … … 1635 1304 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) 1636 1305 { 1637 Info FunctionInfo(__func__); 1638 Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl; 1306 Log() << Verbose(4) << "Adding line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl; 1639 1307 BPS[0] = a; 1640 1308 BPS[1] = b; … … 1644 1312 // increase counter 1645 1313 LinesOnBoundaryCount++; 1646 // also add to open lines1647 CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]);1648 OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT));1649 1314 }; 1650 1315 … … 1654 1319 void Tesselation::AddTesselationTriangle() 1655 1320 { 1656 Info FunctionInfo(__func__);1657 1321 Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1658 1322 … … 1673 1337 void Tesselation::AddTesselationTriangle(const int nr) 1674 1338 { 1675 Info FunctionInfo(__func__); 1676 Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1339 Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1677 1340 1678 1341 // add triangle to global map … … 1692 1355 void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle) 1693 1356 { 1694 Info FunctionInfo(__func__);1695 1357 if (triangle == NULL) 1696 1358 return; 1697 1359 for (int i = 0; i < 3; i++) { 1698 1360 if (triangle->lines[i] != NULL) { 1699 Log() << Verbose( 0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl;1361 Log() << Verbose(5) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl; 1700 1362 triangle->lines[i]->triangles.erase(triangle->Nr); 1701 1363 if (triangle->lines[i]->triangles.empty()) { 1702 Log() << Verbose( 0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl;1364 Log() << Verbose(5) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl; 1703 1365 RemoveTesselationLine(triangle->lines[i]); 1704 1366 } else { 1705 Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: "; 1706 OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL)); 1367 Log() << Verbose(5) << *triangle->lines[i] << " is still attached to another triangle: "; 1707 1368 for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++) 1708 1369 Log() << Verbose(0) << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t"; 1709 1370 Log() << Verbose(0) << endl; 1710 1371 // for (int j=0;j<2;j++) { 1711 // Log() << Verbose( 0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";1372 // Log() << Verbose(5) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": "; 1712 1373 // for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++) 1713 1374 // Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t"; … … 1721 1382 1722 1383 if (TrianglesOnBoundary.erase(triangle->Nr)) 1723 Log() << Verbose( 0) << "Removing triangle Nr. " << triangle->Nr << "." << endl;1384 Log() << Verbose(5) << "Removing triangle Nr. " << triangle->Nr << "." << endl; 1724 1385 delete(triangle); 1725 1386 }; … … 1731 1392 void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line) 1732 1393 { 1733 Info FunctionInfo(__func__);1734 1394 int Numbers[2]; 1735 1395 … … 1752 1412 for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++) 1753 1413 if ((*Runner).second == line) { 1754 Log() << Verbose( 0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;1414 Log() << Verbose(5) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl; 1755 1415 line->endpoints[i]->lines.erase(Runner); 1756 1416 break; … … 1758 1418 } else { // there's just a single line left 1759 1419 if (line->endpoints[i]->lines.erase(line->Nr)) 1760 Log() << Verbose( 0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl;1420 Log() << Verbose(5) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl; 1761 1421 } 1762 1422 if (line->endpoints[i]->lines.empty()) { 1763 Log() << Verbose( 0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl;1423 Log() << Verbose(5) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl; 1764 1424 RemoveTesselationPoint(line->endpoints[i]); 1765 1425 } else { 1766 Log() << Verbose( 0) << *line->endpoints[i] << " has still lines it's attached to: ";1426 Log() << Verbose(5) << *line->endpoints[i] << " has still lines it's attached to: "; 1767 1427 for(LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++) 1768 1428 Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t"; … … 1777 1437 1778 1438 if (LinesOnBoundary.erase(line->Nr)) 1779 Log() << Verbose( 0) << "Removing line Nr. " << line->Nr << "." << endl;1439 Log() << Verbose(5) << "Removing line Nr. " << line->Nr << "." << endl; 1780 1440 delete(line); 1781 1441 }; … … 1788 1448 void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point) 1789 1449 { 1790 Info FunctionInfo(__func__);1791 1450 if (point == NULL) 1792 1451 return; 1793 1452 if (PointsOnBoundary.erase(point->Nr)) 1794 Log() << Verbose( 0) << "Removing point Nr. " << point->Nr << "." << endl;1453 Log() << Verbose(5) << "Removing point Nr. " << point->Nr << "." << endl; 1795 1454 delete(point); 1796 1455 }; … … 1807 1466 int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const 1808 1467 { 1809 Info FunctionInfo(__func__);1810 1468 int adjacentTriangleCount = 0; 1811 1469 class BoundaryPointSet *Points[3]; 1812 1470 1471 Log() << Verbose(2) << "Begin of CheckPresenceOfTriangle" << endl; 1813 1472 // builds a triangle point set (Points) of the end points 1814 1473 for (int i = 0; i < 3; i++) { … … 1829 1488 for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) { 1830 1489 TriangleMap *triangles = &FindLine->second->triangles; 1831 Log() << Verbose( 1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl;1490 Log() << Verbose(3) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl; 1832 1491 for (TriangleMap::const_iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) { 1833 1492 if (FindTriangle->second->IsPresentTupel(Points)) { … … 1835 1494 } 1836 1495 } 1837 Log() << Verbose( 1) << "end." << endl;1496 Log() << Verbose(3) << "end." << endl; 1838 1497 } 1839 1498 // Only one of the triangle lines must be considered for the triangle count. 1840 //Log() << Verbose( 0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;1499 //Log() << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1841 1500 //return adjacentTriangleCount; 1842 1501 } … … 1845 1504 } 1846 1505 1847 Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1506 Log() << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1507 Log() << Verbose(2) << "End of CheckPresenceOfTriangle" << endl; 1848 1508 return adjacentTriangleCount; 1849 1509 }; … … 1859 1519 class BoundaryTriangleSet * Tesselation::GetPresentTriangle(TesselPoint *Candidates[3]) 1860 1520 { 1861 Info FunctionInfo(__func__);1862 1521 class BoundaryTriangleSet *triangle = NULL; 1863 1522 class BoundaryPointSet *Points[3]; … … 1889 1548 } 1890 1549 // Only one of the triangle lines must be considered for the triangle count. 1891 //Log() << Verbose( 0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl;1550 //Log() << Verbose(2) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 1892 1551 //return adjacentTriangleCount; 1893 1552 } … … 1910 1569 void Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC) 1911 1570 { 1912 Info FunctionInfo(__func__);1571 Log() << Verbose(1) << "Begin of FindStartingTriangle\n"; 1913 1572 int i = 0; 1573 TesselPoint* FirstPoint = NULL; 1574 TesselPoint* SecondPoint = NULL; 1914 1575 TesselPoint* MaxPoint[NDIM]; 1915 TesselPoint* Temporary;1916 1576 double maxCoordinate[NDIM]; 1917 BoundaryLineSet BaseLine;1577 Vector Oben; 1918 1578 Vector helper; 1919 1579 Vector Chord; 1920 1580 Vector SearchDirection; 1921 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 1922 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 1923 Vector SphereCenter; 1924 Vector NormalVector; 1925 1926 NormalVector.Zero(); 1581 1582 Oben.Zero(); 1927 1583 1928 1584 for (i = 0; i < 3; i++) { … … 1937 1593 for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) { 1938 1594 const LinkedNodes *List = LC->GetCurrentCell(); 1939 //Log() << Verbose( 1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;1595 //Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 1940 1596 if (List != NULL) { 1941 1597 for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) { 1942 1598 if ((*Runner)->node->x[i] > maxCoordinate[i]) { 1943 Log() << Verbose( 1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl;1599 Log() << Verbose(2) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl; 1944 1600 maxCoordinate[i] = (*Runner)->node->x[i]; 1945 1601 MaxPoint[i] = (*Runner); … … 1952 1608 } 1953 1609 1954 Log() << Verbose( 1) << "Found maximum coordinates: ";1610 Log() << Verbose(2) << "Found maximum coordinates: "; 1955 1611 for (int i=0;i<NDIM;i++) 1956 1612 Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t"; … … 1958 1614 1959 1615 BTS = NULL; 1616 CandidateList *OptCandidates = new CandidateList(); 1960 1617 for (int k=0;k<NDIM;k++) { 1961 NormalVector.Zero();1962 NormalVector.x[k] = 1.;1963 BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]);1964 Log() << Verbose( 0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl;1618 Oben.Zero(); 1619 Oben.x[k] = 1.; 1620 FirstPoint = MaxPoint[k]; 1621 Log() << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl; 1965 1622 1966 1623 double ShortestAngle; 1624 TesselPoint* OptCandidate = NULL; 1967 1625 ShortestAngle = 999999.; // This will contain the angle, which will be always positive (when looking for second point), when looking for third point this will be the quadrant. 1968 1626 1969 FindSecondPointForTesselation(BaseLine.endpoints[0]->node, NormalVector, Temporary, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_... 1970 if (Temporary == NULL) // have we found a second point? 1627 FindSecondPointForTesselation(FirstPoint, Oben, OptCandidate, &ShortestAngle, RADIUS, LC); // we give same point as next candidate as its bonds are looked into in find_second_... 1628 SecondPoint = OptCandidate; 1629 if (SecondPoint == NULL) // have we found a second point? 1971 1630 continue; 1972 BaseLine.endpoints[1] = new BoundaryPointSet(Temporary); 1973 1974 // construct center of circle1975 CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node);1976 CircleCenter.AddVector(BaseLine.endpoints[1]->node->node);1977 CircleCenter.Scale(0.5);1978 1979 // construct normal vector of circle1980 CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node); 1981 C irclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node);1982 1983 double radius = C irclePlaneNormal.NormSquared();1631 1632 helper.CopyVector(FirstPoint->node); 1633 helper.SubtractVector(SecondPoint->node); 1634 helper.Normalize(); 1635 Oben.ProjectOntoPlane(&helper); 1636 Oben.Normalize(); 1637 helper.VectorProduct(&Oben); 1638 ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 1639 1640 Chord.CopyVector(FirstPoint->node); // bring into calling function 1641 Chord.SubtractVector(SecondPoint->node); 1642 double radius = Chord.ScalarProduct(&Chord); 1984 1643 double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.); 1985 1986 NormalVector.ProjectOntoPlane(&CirclePlaneNormal); 1987 NormalVector.Normalize(); 1988 ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 1989 1990 SphereCenter.CopyVector(&NormalVector); 1991 SphereCenter.Scale(CircleRadius); 1992 SphereCenter.AddVector(&CircleCenter); 1993 // Now, NormalVector and SphereCenter are two orthonormalized vectors in the plane defined by CirclePlaneNormal (not normalized) 1644 helper.CopyVector(&Oben); 1645 helper.Scale(CircleRadius); 1646 // Now, oben and helper are two orthonormalized vectors in the plane defined by Chord (not normalized) 1994 1647 1995 1648 // look in one direction of baseline for initial candidate 1996 SearchDirection.MakeNormalVector(&C irclePlaneNormal, &NormalVector); // whether we look "left" first or "right" first is not important ...1649 SearchDirection.MakeNormalVector(&Chord, &Oben); // whether we look "left" first or "right" first is not important ... 1997 1650 1998 1651 // adding point 1 and point 2 and add the line between them 1999 Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl; 2000 Log() << Verbose(0) << "Found second point is at " << *BaseLine.endpoints[1]->node << ".\n"; 2001 2002 //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n"; 2003 CandidateForTesselation OptCandidates(&BaseLine); 2004 FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC); 2005 Log() << Verbose(0) << "List of third Points is:" << endl; 2006 for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) { 2007 Log() << Verbose(0) << " " << *(*it) << endl; 2008 } 2009 2010 BTS = NULL; 2011 AddCandidateTriangle(OptCandidates); 2012 // delete(BaseLine.endpoints[0]); 2013 // delete(BaseLine.endpoints[1]); 2014 1652 Log() << Verbose(1) << "Coordinates of start node at " << *FirstPoint->node << "." << endl; 1653 AddTesselationPoint(FirstPoint, 0); 1654 Log() << Verbose(1) << "Found second point is at " << *SecondPoint->node << ".\n"; 1655 AddTesselationPoint(SecondPoint, 1); 1656 AddTesselationLine(TPS[0], TPS[1], 0); 1657 1658 //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << helper << ".\n"; 1659 FindThirdPointForTesselation(Oben, SearchDirection, helper, BLS[0], NULL, *&OptCandidates, &ShortestAngle, RADIUS, LC); 1660 Log() << Verbose(1) << "List of third Points is "; 1661 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1662 Log() << Verbose(0) << " " << *(*it)->point; 1663 } 1664 Log() << Verbose(0) << endl; 1665 1666 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1667 // add third triangle point 1668 AddTesselationPoint((*it)->point, 2); 1669 // add the second and third line 1670 AddTesselationLine(TPS[1], TPS[2], 1); 1671 AddTesselationLine(TPS[0], TPS[2], 2); 1672 // ... and triangles to the Maps of the Tesselation class 1673 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 1674 AddTesselationTriangle(); 1675 // ... and calculate its normal vector (with correct orientation) 1676 (*it)->OptCenter.Scale(-1.); 1677 Log() << Verbose(2) << "Anti-Oben is currently " << (*it)->OptCenter << "." << endl; 1678 BTS->GetNormalVector((*it)->OptCenter); // vector to compare with should point inwards 1679 Log() << Verbose(0) << "==> Found starting triangle consists of " << *FirstPoint << ", " << *SecondPoint << " and " 1680 << *(*it)->point << " with normal vector " << BTS->NormalVector << ".\n"; 1681 1682 // if we do not reach the end with the next step of iteration, we need to setup a new first line 1683 if (it != OptCandidates->end()--) { 1684 FirstPoint = (*it)->BaseLine->endpoints[0]->node; 1685 SecondPoint = (*it)->point; 1686 // adding point 1 and point 2 and the line between them 1687 AddTesselationPoint(FirstPoint, 0); 1688 AddTesselationPoint(SecondPoint, 1); 1689 AddTesselationLine(TPS[0], TPS[1], 0); 1690 } 1691 Log() << Verbose(2) << "Projection is " << BTS->NormalVector.ScalarProduct(&Oben) << "." << endl; 1692 } 2015 1693 if (BTS != NULL) // we have created one starting triangle 2016 1694 break; 2017 1695 else { 2018 1696 // remove all candidates from the list and then the list itself 2019 OptCandidates.pointlist.clear(); 2020 } 2021 } 1697 class CandidateForTesselation *remover = NULL; 1698 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1699 remover = *it; 1700 delete(remover); 1701 } 1702 OptCandidates->clear(); 1703 } 1704 } 1705 1706 // remove all candidates from the list and then the list itself 1707 class CandidateForTesselation *remover = NULL; 1708 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1709 remover = *it; 1710 delete(remover); 1711 } 1712 delete(OptCandidates); 1713 Log() << Verbose(1) << "End of FindStartingTriangle\n"; 2022 1714 }; 2023 1715 2024 1716 /** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates. 2025 1717 * This is supposed to prevent early closing of the tesselation. 2026 * \param CandidateLine CandidateForTesselation with baseline and shortestangle, i.e. not \a *OptCandidate1718 * \param *BaseRay baseline, i.e. not \a *OptCandidate 2027 1719 * \param *ThirdNode third point in triangle, not in BoundaryLineSet::endpoints 1720 * \param ShortestAngle path length on this circle band for the current \a *ThirdNode 2028 1721 * \param RADIUS radius of sphere 2029 1722 * \param *LC LinkedCell structure 2030 1723 * \return true - there is a better candidate (smaller angle than \a ShortestAngle), false - no better TesselPoint candidate found 2031 1724 */ 2032 //bool Tesselation::HasOtherBaselineBetterCandidate(CandidateForTesselation &CandidateLine, const TesselPoint * const ThirdNode, double RADIUS, const LinkedCell * const LC) const 2033 //{ 2034 // Info FunctionInfo(__func__); 2035 // bool result = false; 2036 // Vector CircleCenter; 2037 // Vector CirclePlaneNormal; 2038 // Vector OldSphereCenter; 2039 // Vector SearchDirection; 2040 // Vector helper; 2041 // TesselPoint *OtherOptCandidate = NULL; 2042 // double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2043 // double radius, CircleRadius; 2044 // BoundaryLineSet *Line = NULL; 2045 // BoundaryTriangleSet *T = NULL; 2046 // 2047 // // check both other lines 2048 // PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr); 2049 // if (FindPoint != PointsOnBoundary.end()) { 2050 // for (int i=0;i<2;i++) { 2051 // LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr); 2052 // if (FindLine != (FindPoint->second)->lines.end()) { 2053 // Line = FindLine->second; 2054 // Log() << Verbose(0) << "Found line " << *Line << "." << endl; 2055 // if (Line->triangles.size() == 1) { 2056 // T = Line->triangles.begin()->second; 2057 // // construct center of circle 2058 // CircleCenter.CopyVector(Line->endpoints[0]->node->node); 2059 // CircleCenter.AddVector(Line->endpoints[1]->node->node); 2060 // CircleCenter.Scale(0.5); 2061 // 2062 // // construct normal vector of circle 2063 // CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node); 2064 // CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node); 2065 // 2066 // // calculate squared radius of circle 2067 // radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2068 // if (radius/4. < RADIUS*RADIUS) { 2069 // CircleRadius = RADIUS*RADIUS - radius/4.; 2070 // CirclePlaneNormal.Normalize(); 2071 // //Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2072 // 2073 // // construct old center 2074 // GetCenterofCircumcircle(&OldSphereCenter, *T->endpoints[0]->node->node, *T->endpoints[1]->node->node, *T->endpoints[2]->node->node); 2075 // helper.CopyVector(&T->NormalVector); // normal vector ensures that this is correct center of the two possible ones 2076 // radius = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter); 2077 // helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2078 // OldSphereCenter.AddVector(&helper); 2079 // OldSphereCenter.SubtractVector(&CircleCenter); 2080 // //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2081 // 2082 // // construct SearchDirection 2083 // SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal); 2084 // helper.CopyVector(Line->endpoints[0]->node->node); 2085 // helper.SubtractVector(ThirdNode->node); 2086 // if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 2087 // SearchDirection.Scale(-1.); 2088 // SearchDirection.ProjectOntoPlane(&OldSphereCenter); 2089 // SearchDirection.Normalize(); 2090 // Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2091 // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 2092 // // rotated the wrong way! 2093 // eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2094 // } 2095 // 2096 // // add third point 2097 // FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, OptCandidates, ThirdNode, RADIUS, LC); 2098 // for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); ++it) { 2099 // if (((*it) == BaseRay->endpoints[0]->node) || ((*it) == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested 2100 // continue; 2101 // Log() << Verbose(0) << " Third point candidate is " << (*it) 2102 // << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 2103 // Log() << Verbose(0) << " Baseline is " << *BaseRay << endl; 2104 // 2105 // // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2106 // TesselPoint *PointCandidates[3]; 2107 // PointCandidates[0] = (*it); 2108 // PointCandidates[1] = BaseRay->endpoints[0]->node; 2109 // PointCandidates[2] = BaseRay->endpoints[1]->node; 2110 // bool check=false; 2111 // int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 2112 // // If there is no triangle, add it regularly. 2113 // if (existentTrianglesCount == 0) { 2114 // SetTesselationPoint((*it), 0); 2115 // SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 2116 // SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 2117 // 2118 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 2119 // OtherOptCandidate = (*it); 2120 // check = true; 2121 // } 2122 // } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 2123 // SetTesselationPoint((*it), 0); 2124 // SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 2125 // SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 2126 // 2127 // // 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) 2128 // // i.e. at least one of the three lines must be present with TriangleCount <= 1 2129 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) { 2130 // OtherOptCandidate = (*it); 2131 // check = true; 2132 // } 2133 // } 2134 // 2135 // if (check) { 2136 // if (ShortestAngle > OtherShortestAngle) { 2137 // Log() << Verbose(0) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl; 2138 // result = true; 2139 // break; 2140 // } 2141 // } 2142 // } 2143 // delete(OptCandidates); 2144 // if (result) 2145 // break; 2146 // } else { 2147 // Log() << Verbose(0) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl; 2148 // } 2149 // } else { 2150 // eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl; 2151 // } 2152 // } else { 2153 // Log() << Verbose(1) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl; 2154 // } 2155 // } 2156 // } else { 2157 // eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl; 2158 // } 2159 // 2160 // return result; 2161 //}; 1725 bool Tesselation::HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const ThirdNode, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const 1726 { 1727 bool result = false; 1728 Vector CircleCenter; 1729 Vector CirclePlaneNormal; 1730 Vector OldSphereCenter; 1731 Vector SearchDirection; 1732 Vector helper; 1733 TesselPoint *OtherOptCandidate = NULL; 1734 double OtherShortestAngle = 2.*M_PI; // This will indicate the quadrant. 1735 double radius, CircleRadius; 1736 BoundaryLineSet *Line = NULL; 1737 BoundaryTriangleSet *T = NULL; 1738 1739 Log() << Verbose(1) << "Begin of HasOtherBaselineBetterCandidate" << endl; 1740 1741 // check both other lines 1742 PointMap::const_iterator FindPoint = PointsOnBoundary.find(ThirdNode->nr); 1743 if (FindPoint != PointsOnBoundary.end()) { 1744 for (int i=0;i<2;i++) { 1745 LineMap::const_iterator FindLine = (FindPoint->second)->lines.find(BaseRay->endpoints[0]->node->nr); 1746 if (FindLine != (FindPoint->second)->lines.end()) { 1747 Line = FindLine->second; 1748 Log() << Verbose(1) << "Found line " << *Line << "." << endl; 1749 if (Line->triangles.size() == 1) { 1750 T = Line->triangles.begin()->second; 1751 // construct center of circle 1752 CircleCenter.CopyVector(Line->endpoints[0]->node->node); 1753 CircleCenter.AddVector(Line->endpoints[1]->node->node); 1754 CircleCenter.Scale(0.5); 1755 1756 // construct normal vector of circle 1757 CirclePlaneNormal.CopyVector(Line->endpoints[0]->node->node); 1758 CirclePlaneNormal.SubtractVector(Line->endpoints[1]->node->node); 1759 1760 // calculate squared radius of circle 1761 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 1762 if (radius/4. < RADIUS*RADIUS) { 1763 CircleRadius = RADIUS*RADIUS - radius/4.; 1764 CirclePlaneNormal.Normalize(); 1765 //Log() << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 1766 1767 // construct old center 1768 GetCenterofCircumcircle(&OldSphereCenter, *T->endpoints[0]->node->node, *T->endpoints[1]->node->node, *T->endpoints[2]->node->node); 1769 helper.CopyVector(&T->NormalVector); // normal vector ensures that this is correct center of the two possible ones 1770 radius = Line->endpoints[0]->node->node->DistanceSquared(&OldSphereCenter); 1771 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 1772 OldSphereCenter.AddVector(&helper); 1773 OldSphereCenter.SubtractVector(&CircleCenter); 1774 //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 1775 1776 // construct SearchDirection 1777 SearchDirection.MakeNormalVector(&T->NormalVector, &CirclePlaneNormal); 1778 helper.CopyVector(Line->endpoints[0]->node->node); 1779 helper.SubtractVector(ThirdNode->node); 1780 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 1781 SearchDirection.Scale(-1.); 1782 SearchDirection.ProjectOntoPlane(&OldSphereCenter); 1783 SearchDirection.Normalize(); 1784 Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 1785 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 1786 // rotated the wrong way! 1787 eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 1788 } 1789 1790 // add third point 1791 CandidateList *OptCandidates = new CandidateList(); 1792 FindThirdPointForTesselation(T->NormalVector, SearchDirection, OldSphereCenter, Line, ThirdNode, OptCandidates, &OtherShortestAngle, RADIUS, LC); 1793 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1794 if (((*it)->point == BaseRay->endpoints[0]->node) || ((*it)->point == BaseRay->endpoints[1]->node)) // skip if it's the same triangle than suggested 1795 continue; 1796 Log() << Verbose(1) << " Third point candidate is " << *(*it)->point 1797 << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 1798 Log() << Verbose(1) << " Baseline is " << *BaseRay << endl; 1799 1800 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 1801 TesselPoint *PointCandidates[3]; 1802 PointCandidates[0] = (*it)->point; 1803 PointCandidates[1] = BaseRay->endpoints[0]->node; 1804 PointCandidates[2] = BaseRay->endpoints[1]->node; 1805 bool check=false; 1806 int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 1807 // If there is no triangle, add it regularly. 1808 if (existentTrianglesCount == 0) { 1809 SetTesselationPoint((*it)->point, 0); 1810 SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 1811 SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 1812 1813 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 1814 OtherOptCandidate = (*it)->point; 1815 check = true; 1816 } 1817 } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 1818 SetTesselationPoint((*it)->point, 0); 1819 SetTesselationPoint(BaseRay->endpoints[0]->node, 1); 1820 SetTesselationPoint(BaseRay->endpoints[1]->node, 2); 1821 1822 // 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) 1823 // i.e. at least one of the three lines must be present with TriangleCount <= 1 1824 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) { 1825 OtherOptCandidate = (*it)->point; 1826 check = true; 1827 } 1828 } 1829 1830 if (check) { 1831 if (ShortestAngle > OtherShortestAngle) { 1832 Log() << Verbose(1) << "There is a better candidate than " << *ThirdNode << " with " << ShortestAngle << " from baseline " << *Line << ": " << *OtherOptCandidate << " with " << OtherShortestAngle << "." << endl; 1833 result = true; 1834 break; 1835 } 1836 } 1837 } 1838 delete(OptCandidates); 1839 if (result) 1840 break; 1841 } else { 1842 Log() << Verbose(1) << "Circumcircle for base line " << *Line << " and base triangle " << T << " is too big!" << endl; 1843 } 1844 } else { 1845 eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl; 1846 } 1847 } else { 1848 Log() << Verbose(2) << "No present baseline between " << BaseRay->endpoints[0] << " and candidate " << *ThirdNode << "." << endl; 1849 } 1850 } 1851 } else { 1852 eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl; 1853 } 1854 1855 Log() << Verbose(1) << "End of HasOtherBaselineBetterCandidate" << endl; 1856 1857 return result; 1858 }; 2162 1859 2163 1860 /** This function finds a triangle to a line, adjacent to an existing one. 2164 1861 * @param out output stream for debugging 2165 * @param CandidateLine current cadndiatebaseline to search from1862 * @param Line current baseline to search from 2166 1863 * @param T current triangle which \a Line is edge of 2167 1864 * @param RADIUS radius of the rolling ball … … 2169 1866 * @param *LC LinkedCell structure with neighbouring points 2170 1867 */ 2171 bool Tesselation::FindNextSuitableTriangle( CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC)2172 { 2173 Info FunctionInfo(__func__);1868 bool Tesselation::FindNextSuitableTriangle(BoundaryLineSet &Line, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC) 1869 { 1870 Log() << Verbose(0) << "Begin of FindNextSuitableTriangle\n"; 2174 1871 bool result = true; 1872 CandidateList *OptCandidates = new CandidateList(); 2175 1873 2176 1874 Vector CircleCenter; 2177 1875 Vector CirclePlaneNormal; 2178 Vector RelativeSphereCenter;1876 Vector OldSphereCenter; 2179 1877 Vector SearchDirection; 2180 1878 Vector helper; 2181 1879 TesselPoint *ThirdNode = NULL; 2182 1880 LineMap::iterator testline; 1881 double ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2183 1882 double radius, CircleRadius; 2184 1883 1884 Log() << Verbose(1) << "Current baseline is " << Line << " of triangle " << T << "." << endl; 2185 1885 for (int i=0;i<3;i++) 2186 if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) {1886 if ((T.endpoints[i]->node != Line.endpoints[0]->node) && (T.endpoints[i]->node != Line.endpoints[1]->node)) 2187 1887 ThirdNode = T.endpoints[i]->node; 2188 break;2189 }2190 Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdNode " << *ThirdNode << " of triangle " << T << "." << endl;2191 1888 2192 1889 // construct center of circle 2193 CircleCenter.CopyVector( CandidateLine.BaseLine->endpoints[0]->node->node);2194 CircleCenter.AddVector( CandidateLine.BaseLine->endpoints[1]->node->node);1890 CircleCenter.CopyVector(Line.endpoints[0]->node->node); 1891 CircleCenter.AddVector(Line.endpoints[1]->node->node); 2195 1892 CircleCenter.Scale(0.5); 2196 1893 2197 1894 // construct normal vector of circle 2198 CirclePlaneNormal.CopyVector( CandidateLine.BaseLine->endpoints[0]->node->node);2199 CirclePlaneNormal.SubtractVector( CandidateLine.BaseLine->endpoints[1]->node->node);1895 CirclePlaneNormal.CopyVector(Line.endpoints[0]->node->node); 1896 CirclePlaneNormal.SubtractVector(Line.endpoints[1]->node->node); 2200 1897 2201 1898 // calculate squared radius of circle 2202 1899 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2203 1900 if (radius/4. < RADIUS*RADIUS) { 2204 // construct relative sphere center with now known CircleCenter2205 RelativeSphereCenter.CopyVector(&T.SphereCenter);2206 RelativeSphereCenter.SubtractVector(&CircleCenter);2207 2208 1901 CircleRadius = RADIUS*RADIUS - radius/4.; 2209 1902 CirclePlaneNormal.Normalize(); 2210 Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2211 2212 Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl; 2213 2214 // construct SearchDirection and an "outward pointer" 2215 SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal); 2216 helper.CopyVector(&CircleCenter); 1903 //Log() << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 1904 1905 // construct old center 1906 GetCenterofCircumcircle(&OldSphereCenter, *T.endpoints[0]->node->node, *T.endpoints[1]->node->node, *T.endpoints[2]->node->node); 1907 helper.CopyVector(&T.NormalVector); // normal vector ensures that this is correct center of the two possible ones 1908 radius = Line.endpoints[0]->node->node->DistanceSquared(&OldSphereCenter); 1909 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 1910 OldSphereCenter.AddVector(&helper); 1911 OldSphereCenter.SubtractVector(&CircleCenter); 1912 //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 1913 1914 // construct SearchDirection 1915 SearchDirection.MakeNormalVector(&T.NormalVector, &CirclePlaneNormal); 1916 helper.CopyVector(Line.endpoints[0]->node->node); 2217 1917 helper.SubtractVector(ThirdNode->node); 2218 1918 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 2219 1919 SearchDirection.Scale(-1.); 2220 Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2221 if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 1920 SearchDirection.ProjectOntoPlane(&OldSphereCenter); 1921 SearchDirection.Normalize(); 1922 Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 1923 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 2222 1924 // rotated the wrong way! 2223 1925 eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; … … 2225 1927 2226 1928 // add third point 2227 FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC);1929 FindThirdPointForTesselation(T.NormalVector, SearchDirection, OldSphereCenter, &Line, ThirdNode, OptCandidates, &ShortestAngle, RADIUS, LC); 2228 1930 2229 1931 } else { 2230 Log() << Verbose( 0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl;2231 } 2232 2233 if ( CandidateLine.pointlist.empty()) {1932 Log() << Verbose(1) << "Circumcircle for base line " << Line << " and base triangle " << T << " is too big!" << endl; 1933 } 1934 1935 if (OptCandidates->begin() == OptCandidates->end()) { 2234 1936 eLog() << Verbose(2) << "Could not find a suitable candidate." << endl; 2235 1937 return false; 2236 1938 } 2237 Log() << Verbose(0) << "Third Points are: " << endl; 2238 for (TesselPointList::iterator it = CandidateLine.pointlist.begin(); it != CandidateLine.pointlist.end(); ++it) { 2239 Log() << Verbose(0) << " " << *(*it) << endl; 2240 } 2241 2242 return true; 2243 2244 // BoundaryLineSet *BaseRay = CandidateLine.BaseLine; 2245 // for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2246 // Log() << Verbose(0) << "Third point candidate is " << *(*it)->point 2247 // << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 2248 // Log() << Verbose(0) << "Baseline is " << *BaseRay << endl; 2249 // 2250 // // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2251 // TesselPoint *PointCandidates[3]; 2252 // PointCandidates[0] = (*it)->point; 2253 // PointCandidates[1] = BaseRay->endpoints[0]->node; 2254 // PointCandidates[2] = BaseRay->endpoints[1]->node; 2255 // int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 2256 // 2257 // BTS = NULL; 2258 // // check for present edges and whether we reach better candidates from them 2259 // //if (HasOtherBaselineBetterCandidate(BaseRay, (*it)->point, ShortestAngle, RADIUS, LC) ) { 2260 // if (0) { 2261 // result = false; 2262 // break; 2263 // } else { 2264 // // If there is no triangle, add it regularly. 2265 // if (existentTrianglesCount == 0) { 2266 // AddTesselationPoint((*it)->point, 0); 2267 // AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2268 // AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2269 // 2270 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 2271 // CandidateLine.point = (*it)->point; 2272 // CandidateLine.OptCenter.CopyVector(&((*it)->OptCenter)); 2273 // CandidateLine.OtherOptCenter.CopyVector(&((*it)->OtherOptCenter)); 2274 // CandidateLine.ShortestAngle = ShortestAngle; 2275 // } else { 2276 //// eLog() << Verbose(1) << "This triangle consisting of "; 2277 //// Log() << Verbose(0) << *(*it)->point << ", "; 2278 //// Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2279 //// Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2280 //// Log() << Verbose(0) << "exists and is not added, as it 0x80000000006fc150(does not seem helpful!" << endl; 2281 // result = false; 2282 // } 2283 // } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 2284 // AddTesselationPoint((*it)->point, 0); 2285 // AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2286 // AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2287 // 2288 // // 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) 2289 // // i.e. at least one of the three lines must be present with TriangleCount <= 1 2290 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS) || CandidateLine.BaseLine->skipped) { 2291 // CandidateLine.point = (*it)->point; 2292 // CandidateLine.OptCenter.CopyVector(&(*it)->OptCenter); 2293 // CandidateLine.OtherOptCenter.CopyVector(&(*it)->OtherOptCenter); 2294 // CandidateLine.ShortestAngle = ShortestAngle+2.*M_PI; 2295 // 2296 // } else { 2297 //// 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; 2298 // result = false; 2299 // } 2300 // } else { 2301 //// Log() << Verbose(1) << "This triangle consisting of "; 2302 //// Log() << Verbose(0) << *(*it)->point << ", "; 2303 //// Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2304 //// Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2305 //// Log() << Verbose(0) << "is invalid!" << endl; 2306 // result = false; 2307 // } 2308 // } 2309 // 2310 // // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 2311 // BaseRay = BLS[0]; 2312 // if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) { 2313 // eLog() << Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl; 2314 // exit(255); 2315 // } 2316 // 2317 // } 2318 // 2319 // // remove all candidates from the list and then the list itself 2320 // class CandidateForTesselation *remover = NULL; 2321 // for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2322 // remover = *it; 2323 // delete(remover); 2324 // } 2325 // delete(OptCandidates); 1939 Log() << Verbose(1) << "Third Points are "; 1940 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1941 Log() << Verbose(1) << " " << *(*it)->point << endl; 1942 } 1943 1944 BoundaryLineSet *BaseRay = &Line; 1945 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 1946 Log() << Verbose(1) << " Third point candidate is " << *(*it)->point 1947 << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 1948 Log() << Verbose(1) << " Baseline is " << *BaseRay << endl; 1949 1950 // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 1951 TesselPoint *PointCandidates[3]; 1952 PointCandidates[0] = (*it)->point; 1953 PointCandidates[1] = BaseRay->endpoints[0]->node; 1954 PointCandidates[2] = BaseRay->endpoints[1]->node; 1955 int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 1956 1957 BTS = NULL; 1958 // check for present edges and whether we reach better candidates from them 1959 if (HasOtherBaselineBetterCandidate(BaseRay, (*it)->point, ShortestAngle, RADIUS, LC) ) { 1960 result = false; 1961 break; 1962 } else { 1963 // If there is no triangle, add it regularly. 1964 if (existentTrianglesCount == 0) { 1965 AddTesselationPoint((*it)->point, 0); 1966 AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 1967 AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 1968 1969 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 1970 AddTesselationLine(TPS[0], TPS[1], 0); 1971 AddTesselationLine(TPS[0], TPS[2], 1); 1972 AddTesselationLine(TPS[1], TPS[2], 2); 1973 1974 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 1975 AddTesselationTriangle(); 1976 (*it)->OptCenter.Scale(-1.); 1977 BTS->GetNormalVector((*it)->OptCenter); 1978 (*it)->OptCenter.Scale(-1.); 1979 1980 Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector 1981 << " for this triangle ... " << endl; 1982 //Log() << Verbose(1) << "We have "<< TrianglesOnBoundaryCount << " for line " << *BaseRay << "." << endl; 1983 } else { 1984 eLog() << Verbose(2) << "This triangle consisting of "; 1985 Log() << Verbose(0) << *(*it)->point << ", "; 1986 Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 1987 Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 1988 Log() << Verbose(0) << "exists and is not added, as it does not seem helpful!" << endl; 1989 result = false; 1990 } 1991 } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 1992 AddTesselationPoint((*it)->point, 0); 1993 AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 1994 AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 1995 1996 // 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) 1997 // i.e. at least one of the three lines must be present with TriangleCount <= 1 1998 if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS)) { 1999 AddTesselationLine(TPS[0], TPS[1], 0); 2000 AddTesselationLine(TPS[0], TPS[2], 1); 2001 AddTesselationLine(TPS[1], TPS[2], 2); 2002 2003 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2004 AddTesselationTriangle(); // add to global map 2005 2006 (*it)->OtherOptCenter.Scale(-1.); 2007 BTS->GetNormalVector((*it)->OtherOptCenter); 2008 (*it)->OtherOptCenter.Scale(-1.); 2009 2010 eLog() << Verbose(2) << "--> WARNING: Special new triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " for this triangle ... " << endl; 2011 Log() << Verbose(1) << "We have "<< BaseRay->triangles.size() << " for line " << BaseRay << "." << endl; 2012 } else { 2013 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; 2014 result = false; 2015 } 2016 } else { 2017 Log() << Verbose(1) << "This triangle consisting of "; 2018 Log() << Verbose(0) << *(*it)->point << ", "; 2019 Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2020 Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2021 Log() << Verbose(0) << "is invalid!" << endl; 2022 result = false; 2023 } 2024 } 2025 2026 // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 2027 BaseRay = BLS[0]; 2028 if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) { 2029 eLog() << Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl; 2030 exit(255); 2031 } 2032 2033 } 2034 2035 // remove all candidates from the list and then the list itself 2036 class CandidateForTesselation *remover = NULL; 2037 for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2038 remover = *it; 2039 delete(remover); 2040 } 2041 delete(OptCandidates); 2042 Log() << Verbose(0) << "End of FindNextSuitableTriangle\n"; 2326 2043 return result; 2327 };2328 2329 /** Adds the present line and candidate point from \a &CandidateLine to the Tesselation.2330 * \param CandidateLine triangle to add2331 * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in AddTesselationLine()2332 */2333 void Tesselation::AddCandidateTriangle(CandidateForTesselation CandidateLine)2334 {2335 Info FunctionInfo(__func__);2336 Vector Center;2337 TesselPoint * const TurningPoint = CandidateLine.BaseLine->endpoints[0]->node;2338 2339 // fill the set of neighbours2340 Center.CopyVector(CandidateLine.BaseLine->endpoints[1]->node->node);2341 Center.SubtractVector(TurningPoint->node);2342 set<TesselPoint*> SetOfNeighbours;2343 SetOfNeighbours.insert(CandidateLine.BaseLine->endpoints[1]->node);2344 for (TesselPointList::iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); Runner++)2345 SetOfNeighbours.insert(*Runner);2346 TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, &Center);2347 2348 // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles)2349 TesselPointList::iterator Runner = connectedClosestPoints->begin();2350 TesselPointList::iterator Sprinter = Runner;2351 Sprinter++;2352 while(Sprinter != connectedClosestPoints->end()) {2353 // add the points2354 AddTesselationPoint(TurningPoint, 0);2355 AddTesselationPoint((*Runner), 1);2356 AddTesselationPoint((*Sprinter), 2);2357 2358 2359 // add the lines2360 AddTesselationLine(TPS[0], TPS[1], 0);2361 AddTesselationLine(TPS[0], TPS[2], 1);2362 AddTesselationLine(TPS[1], TPS[2], 2);2363 2364 // add the triangles2365 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount);2366 AddTesselationTriangle();2367 BTS->GetCenter(&Center);2368 Center.SubtractVector(&CandidateLine.OptCenter);2369 BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter);2370 BTS->GetNormalVector(Center);2371 2372 Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << "." << endl;2373 Runner = Sprinter;2374 Sprinter++;2375 }2376 delete(connectedClosestPoints);2377 2044 }; 2378 2045 … … 2386 2053 class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base) 2387 2054 { 2388 Info FunctionInfo(__func__);2389 2055 class BoundaryPointSet *Spot = NULL; 2390 2056 class BoundaryLineSet *OtherBase; … … 2398 2064 OtherBase = new class BoundaryLineSet(BPS,-1); 2399 2065 2400 Log() << Verbose( 1) << "INFO: Current base line is " << *Base << "." << endl;2401 Log() << Verbose( 1) << "INFO: Other base line is " << *OtherBase << "." << endl;2066 Log() << Verbose(3) << "INFO: Current base line is " << *Base << "." << endl; 2067 Log() << Verbose(3) << "INFO: Other base line is " << *OtherBase << "." << endl; 2402 2068 2403 2069 // get the closest point on each line to the other line … … 2419 2085 delete(ClosestPoint); 2420 2086 if ((distance[0] * distance[1]) > 0) { // have same sign? 2421 Log() << Verbose( 1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl;2087 Log() << Verbose(3) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl; 2422 2088 if (distance[0] < distance[1]) { 2423 2089 Spot = Base->endpoints[0]; … … 2427 2093 return Spot; 2428 2094 } else { // different sign, i.e. we are in between 2429 Log() << Verbose( 0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl;2095 Log() << Verbose(3) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl; 2430 2096 return NULL; 2431 2097 } … … 2435 2101 void Tesselation::PrintAllBoundaryPoints(ofstream *out) const 2436 2102 { 2437 Info FunctionInfo(__func__);2438 2103 // print all lines 2439 Log() << Verbose( 0) << "Printing all boundary points for debugging:" << endl;2104 Log() << Verbose(1) << "Printing all boundary points for debugging:" << endl; 2440 2105 for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++) 2441 Log() << Verbose( 0) << *(PointRunner->second) << endl;2106 Log() << Verbose(2) << *(PointRunner->second) << endl; 2442 2107 }; 2443 2108 2444 2109 void Tesselation::PrintAllBoundaryLines(ofstream *out) const 2445 2110 { 2446 Info FunctionInfo(__func__);2447 2111 // print all lines 2448 Log() << Verbose( 0) << "Printing all boundary lines for debugging:" << endl;2112 Log() << Verbose(1) << "Printing all boundary lines for debugging:" << endl; 2449 2113 for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++) 2450 Log() << Verbose( 0) << *(LineRunner->second) << endl;2114 Log() << Verbose(2) << *(LineRunner->second) << endl; 2451 2115 }; 2452 2116 2453 2117 void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const 2454 2118 { 2455 Info FunctionInfo(__func__);2456 2119 // print all triangles 2457 Log() << Verbose( 0) << "Printing all boundary triangles for debugging:" << endl;2120 Log() << Verbose(1) << "Printing all boundary triangles for debugging:" << endl; 2458 2121 for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++) 2459 Log() << Verbose( 0) << *(TriangleRunner->second) << endl;2122 Log() << Verbose(2) << *(TriangleRunner->second) << endl; 2460 2123 }; 2461 2124 … … 2467 2130 double Tesselation::PickFarthestofTwoBaselines(class BoundaryLineSet *Base) 2468 2131 { 2469 Info FunctionInfo(__func__);2470 2132 class BoundaryLineSet *OtherBase; 2471 2133 Vector *ClosestPoint[2]; … … 2479 2141 OtherBase = new class BoundaryLineSet(BPS,-1); 2480 2142 2481 Log() << Verbose( 0) << "INFO: Current base line is " << *Base << "." << endl;2482 Log() << Verbose( 0) << "INFO: Other base line is " << *OtherBase << "." << endl;2143 Log() << Verbose(3) << "INFO: Current base line is " << *Base << "." << endl; 2144 Log() << Verbose(3) << "INFO: Other base line is " << *OtherBase << "." << endl; 2483 2145 2484 2146 // get the closest point on each line to the other line … … 2500 2162 2501 2163 if (Distance.NormSquared() < MYEPSILON) { // check for intersection 2502 Log() << Verbose( 0) << "REJECT: Both lines have an intersection: Nothing to do." << endl;2164 Log() << Verbose(3) << "REJECT: Both lines have an intersection: Nothing to do." << endl; 2503 2165 return false; 2504 2166 } else { // check for sign against BaseLineNormal … … 2510 2172 } 2511 2173 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2512 Log() << Verbose( 1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;2174 Log() << Verbose(4) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl; 2513 2175 BaseLineNormal.AddVector(&(runner->second->NormalVector)); 2514 2176 } … … 2516 2178 2517 2179 if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip 2518 Log() << Verbose( 0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl;2180 Log() << Verbose(2) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl; 2519 2181 // calculate volume summand as a general tetraeder 2520 2182 return volume; 2521 2183 } else { // Base higher than OtherBase -> do nothing 2522 Log() << Verbose( 0) << "REJECT: Base line is higher: Nothing to do." << endl;2184 Log() << Verbose(2) << "REJECT: Base line is higher: Nothing to do." << endl; 2523 2185 return 0.; 2524 2186 } … … 2535 2197 class BoundaryLineSet * Tesselation::FlipBaseline(class BoundaryLineSet *Base) 2536 2198 { 2537 Info FunctionInfo(__func__);2538 2199 class BoundaryLineSet *OldLines[4], *NewLine; 2539 2200 class BoundaryPointSet *OldPoints[2]; … … 2542 2203 int i,m; 2543 2204 2205 Log() << Verbose(1) << "Begin of FlipBaseline" << endl; 2206 2544 2207 // calculate NormalVector for later use 2545 2208 BaseLineNormal.Zero(); … … 2549 2212 } 2550 2213 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2551 Log() << Verbose( 1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl;2214 Log() << Verbose(4) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl; 2552 2215 BaseLineNormal.AddVector(&(runner->second->NormalVector)); 2553 2216 } … … 2562 2225 i=0; 2563 2226 m=0; 2564 Log() << Verbose( 0) << "The four old lines are: ";2227 Log() << Verbose(3) << "The four old lines are: "; 2565 2228 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2566 2229 for (int j=0;j<3;j++) // all of their endpoints and baselines … … 2570 2233 } 2571 2234 Log() << Verbose(0) << endl; 2572 Log() << Verbose( 0) << "The two old points are: ";2235 Log() << Verbose(3) << "The two old points are: "; 2573 2236 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2574 2237 for (int j=0;j<3;j++) // all of their endpoints and baselines … … 2596 2259 2597 2260 // remove triangles and baseline removes itself 2598 Log() << Verbose( 0) << "INFO: Deleting baseline " << *Base << " from global list." << endl;2261 Log() << Verbose(3) << "INFO: Deleting baseline " << *Base << " from global list." << endl; 2599 2262 OldBaseLineNr = Base->Nr; 2600 2263 m=0; 2601 2264 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2602 Log() << Verbose( 0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl;2265 Log() << Verbose(3) << "INFO: Deleting triangle " << *(runner->second) << "." << endl; 2603 2266 OldTriangleNrs[m++] = runner->second->Nr; 2604 2267 RemoveTesselationTriangle(runner->second); … … 2610 2273 NewLine = new class BoundaryLineSet(BPS, OldBaseLineNr); 2611 2274 LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one 2612 Log() << Verbose( 0) << "INFO: Created new baseline " << *NewLine << "." << endl;2275 Log() << Verbose(3) << "INFO: Created new baseline " << *NewLine << "." << endl; 2613 2276 2614 2277 // construct new triangles with flipped baseline … … 2625 2288 BTS->GetNormalVector(BaseLineNormal); 2626 2289 AddTesselationTriangle(OldTriangleNrs[0]); 2627 Log() << Verbose( 0) << "INFO: Created new triangle " << *BTS << "." << endl;2290 Log() << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl; 2628 2291 2629 2292 BLS[0] = (i==2 ? OldLines[3] : OldLines[2]); … … 2633 2296 BTS->GetNormalVector(BaseLineNormal); 2634 2297 AddTesselationTriangle(OldTriangleNrs[1]); 2635 Log() << Verbose( 0) << "INFO: Created new triangle " << *BTS << "." << endl;2298 Log() << Verbose(3) << "INFO: Created new triangle " << *BTS << "." << endl; 2636 2299 } else { 2637 eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl;2300 Log() << Verbose(1) << "The four old lines do not connect, something's utterly wrong here!" << endl; 2638 2301 return NULL; 2639 2302 } 2640 2303 2304 Log() << Verbose(1) << "End of FlipBaseline" << endl; 2641 2305 return NewLine; 2642 2306 }; … … 2653 2317 void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC) 2654 2318 { 2655 Info FunctionInfo(__func__);2319 Log() << Verbose(2) << "Begin of FindSecondPointForTesselation" << endl; 2656 2320 Vector AngleCheck; 2657 2321 class TesselPoint* Candidate = NULL; … … 2674 2338 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2675 2339 } 2676 Log() << Verbose( 0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :"2340 Log() << Verbose(3) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" 2677 2341 << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl; 2678 2342 … … 2681 2345 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2682 2346 const LinkedNodes *List = LC->GetCurrentCell(); 2683 //Log() << Verbose( 1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;2347 //Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2684 2348 if (List != NULL) { 2685 2349 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 2712 2376 angle = AngleCheck.Angle(&Oben); 2713 2377 if (angle < Storage[0]) { 2714 //Log() << Verbose( 1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]);2715 Log() << Verbose( 1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n";2378 //Log() << Verbose(3) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]); 2379 Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n"; 2716 2380 OptCandidate = Candidate; 2717 2381 Storage[0] = angle; 2718 //Log() << Verbose( 1) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]);2382 //Log() << Verbose(3) << "Changing something in Storage: %lf %lf. \n", Storage[0], Storage[2]); 2719 2383 } else { 2720 //Log() << Verbose( 1) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl;2384 //Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Looses with angle " << angle << " to a better candidate " << *OptCandidate << endl; 2721 2385 } 2722 2386 } else { 2723 //Log() << Verbose( 1) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl;2387 //Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Refused due to Radius " << norm << endl; 2724 2388 } 2725 2389 } else { 2726 //Log() << Verbose( 1) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl;2390 //Log() << Verbose(3) << "Current candidate is " << *Candidate << ": Candidate is equal to first endpoint." << *a << "." << endl; 2727 2391 } 2728 2392 } 2729 2393 } else { 2730 Log() << Verbose( 0) << "Linked cell list is empty." << endl;2394 Log() << Verbose(3) << "Linked cell list is empty." << endl; 2731 2395 } 2732 2396 } 2397 Log() << Verbose(2) << "End of FindSecondPointForTesselation" << endl; 2733 2398 }; 2734 2399 … … 2759 2424 * @param SearchDirection general direction where to search for the next point, relative to center of BaseLine 2760 2425 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle 2761 * @param CandidateLine CandidateForTesselation with the current base line and list of candidates and ShortestAngle2426 * @param BaseLine BoundaryLineSet with the current base line 2762 2427 * @param ThirdNode third point to avoid in search 2428 * @param candidates list of equally good candidates to return 2429 * @param ShortestAngle the current path length on this circle band for the current OptCandidate 2763 2430 * @param RADIUS radius of sphere 2764 2431 * @param *LC LinkedCell structure with neighbouring points 2765 2432 */ 2766 void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class TesselPoint * const ThirdNode, const double RADIUS, const LinkedCell *LC) const 2767 { 2768 Info FunctionInfo(__func__); 2433 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 2434 { 2769 2435 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 2770 2436 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in … … 2774 2440 Vector NewNormalVector; // normal vector of the Candidate's triangle 2775 2441 Vector helper, OptCandidateCenter, OtherOptCandidateCenter; 2776 Vector RelativeOldSphereCenter;2777 Vector NewPlaneCenter;2778 2442 double CircleRadius; // radius of this circle 2779 2443 double radius; 2780 double otherradius;2781 2444 double alpha, Otheralpha; // angles (i.e. parameter for the circle). 2782 2445 int N[NDIM], Nlower[NDIM], Nupper[NDIM]; 2783 2446 TesselPoint *Candidate = NULL; 2784 2785 Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 2447 CandidateForTesselation *optCandidate = NULL; 2448 2449 Log() << Verbose(1) << "Begin of FindThirdPointForTesselation" << endl; 2450 2451 Log() << Verbose(2) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 2786 2452 2787 2453 // construct center of circle 2788 CircleCenter.CopyVector( CandidateLine.BaseLine->endpoints[0]->node->node);2789 CircleCenter.AddVector( CandidateLine.BaseLine->endpoints[1]->node->node);2454 CircleCenter.CopyVector(BaseLine->endpoints[0]->node->node); 2455 CircleCenter.AddVector(BaseLine->endpoints[1]->node->node); 2790 2456 CircleCenter.Scale(0.5); 2791 2457 2792 2458 // construct normal vector of circle 2793 CirclePlaneNormal.CopyVector(CandidateLine.BaseLine->endpoints[0]->node->node); 2794 CirclePlaneNormal.SubtractVector(CandidateLine.BaseLine->endpoints[1]->node->node); 2795 2796 RelativeOldSphereCenter.CopyVector(&OldSphereCenter); 2797 RelativeOldSphereCenter.SubtractVector(&CircleCenter); 2459 CirclePlaneNormal.CopyVector(BaseLine->endpoints[0]->node->node); 2460 CirclePlaneNormal.SubtractVector(BaseLine->endpoints[1]->node->node); 2798 2461 2799 2462 // calculate squared radius TesselPoint *ThirdNode,f circle 2800 radius = CirclePlaneNormal. NormSquared()/4.;2801 if (radius < RADIUS*RADIUS) {2802 CircleRadius = RADIUS*RADIUS - radius ;2463 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2464 if (radius/4. < RADIUS*RADIUS) { 2465 CircleRadius = RADIUS*RADIUS - radius/4.; 2803 2466 CirclePlaneNormal.Normalize(); 2804 Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl;2467 //Log() << Verbose(2) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2805 2468 2806 2469 // test whether old center is on the band's plane 2807 if (fabs( RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {2808 eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;2809 RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal);2810 } 2811 radius = RelativeOldSphereCenter.NormSquared();2470 if (fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 2471 eLog() << Verbose(1) << "Something's very wrong here: OldSphereCenter is not on the band's plane as desired by " << fabs(OldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 2472 OldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal); 2473 } 2474 radius = OldSphereCenter.ScalarProduct(&OldSphereCenter); 2812 2475 if (fabs(radius - CircleRadius) < HULLEPSILON) { 2813 Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl;2476 //Log() << Verbose(2) << "INFO: OldSphereCenter is at " << OldSphereCenter << "." << endl; 2814 2477 2815 2478 // check SearchDirection 2816 Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl;2817 if (fabs( RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!2479 //Log() << Verbose(2) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2480 if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2818 2481 eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl; 2819 2482 } … … 2823 2486 for(int i=0;i<NDIM;i++) // store indices of this cell 2824 2487 N[i] = LC->n[i]; 2825 //Log() << Verbose( 1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;2488 //Log() << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 2826 2489 } else { 2827 2490 eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; … … 2829 2492 } 2830 2493 // then go through the current and all neighbouring cells and check the contained points for possible candidates 2831 //Log() << Verbose( 1) << "LC Intervals:";2494 //Log() << Verbose(2) << "LC Intervals:"; 2832 2495 for (int i=0;i<NDIM;i++) { 2833 2496 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; … … 2840 2503 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 2841 2504 const LinkedNodes *List = LC->GetCurrentCell(); 2842 //Log() << Verbose( 1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl;2505 //Log() << Verbose(2) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2843 2506 if (List != NULL) { 2844 2507 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 2846 2509 2847 2510 // check for three unique points 2848 Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter<< "." << endl;2849 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){2850 2851 // find center on the plane2852 GetCenterofCircumcircle(&New PlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node);2853 Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl;2854 2855 if ( NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node)2856 && (fabs(NewNormalVector. NormSquared()) > HULLEPSILON)2511 //Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " at " << Candidate->node << "." << endl; 2512 if ((Candidate != BaseLine->endpoints[0]->node) && (Candidate != BaseLine->endpoints[1]->node) ){ 2513 2514 // construct both new centers 2515 GetCenterofCircumcircle(&NewSphereCenter, *BaseLine->endpoints[0]->node->node, *BaseLine->endpoints[1]->node->node, *Candidate->node); 2516 OtherNewSphereCenter.CopyVector(&NewSphereCenter); 2517 2518 if ((NewNormalVector.MakeNormalVector(BaseLine->endpoints[0]->node->node, BaseLine->endpoints[1]->node->node, Candidate->node)) 2519 && (fabs(NewNormalVector.ScalarProduct(&NewNormalVector)) > HULLEPSILON) 2857 2520 ) { 2858 Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2859 radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter); 2860 Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2861 Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2862 Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl; 2521 helper.CopyVector(&NewNormalVector); 2522 //Log() << Verbose(2) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 2523 radius = BaseLine->endpoints[0]->node->node->DistanceSquared(&NewSphereCenter); 2863 2524 if (radius < RADIUS*RADIUS) { 2864 otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter);2865 if (fabs(radius - otherradius) > HULLEPSILON) {2866 eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl;2867 }2868 // construct both new centers2869 NewSphereCenter.CopyVector(&NewPlaneCenter);2870 OtherNewSphereCenter.CopyVector(&NewPlaneCenter);2871 helper.CopyVector(&NewNormalVector);2872 2525 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2873 Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper<< " with sphere radius " << RADIUS << "." << endl;2526 //Log() << Verbose(2) << "INFO: Distance of NewCircleCenter to NewSphereCenter is " << helper.Norm() << " with sphere radius " << RADIUS << "." << endl; 2874 2527 NewSphereCenter.AddVector(&helper); 2875 Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2528 NewSphereCenter.SubtractVector(&CircleCenter); 2529 //Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 2530 2876 2531 // OtherNewSphereCenter is created by the same vector just in the other direction 2877 2532 helper.Scale(-1.); 2878 2533 OtherNewSphereCenter.AddVector(&helper); 2879 Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2534 OtherNewSphereCenter.SubtractVector(&CircleCenter); 2535 //Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 2880 2536 2881 2537 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2882 2538 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 2883 2539 alpha = min(alpha, Otheralpha); 2884 2885 2540 // if there is a better candidate, drop the current list and add the new candidate 2886 2541 // otherwise ignore the new candidate and keep the list 2887 if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) { 2542 if (*ShortestAngle > (alpha - HULLEPSILON)) { 2543 optCandidate = new CandidateForTesselation(Candidate, BaseLine, OptCandidateCenter, OtherOptCandidateCenter); 2888 2544 if (fabs(alpha - Otheralpha) > MYEPSILON) { 2889 CandidateLine.OptCenter.CopyVector(&NewSphereCenter);2890 CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter);2545 optCandidate->OptCenter.CopyVector(&NewSphereCenter); 2546 optCandidate->OtherOptCenter.CopyVector(&OtherNewSphereCenter); 2891 2547 } else { 2892 CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter);2893 CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter);2548 optCandidate->OptCenter.CopyVector(&OtherNewSphereCenter); 2549 optCandidate->OtherOptCenter.CopyVector(&NewSphereCenter); 2894 2550 } 2895 2551 // if there is an equal candidate, add it to the list without clearing the list 2896 if (( CandidateLine.ShortestAngle - HULLEPSILON) < alpha) {2897 CandidateLine.pointlist.push_back(Candidate);2898 Log() << Verbose( 0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with "2899 << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl;2552 if ((*ShortestAngle - HULLEPSILON) < alpha) { 2553 candidates->push_back(optCandidate); 2554 Log() << Verbose(2) << "ACCEPT: We have found an equally good candidate: " << *(optCandidate->point) << " with " 2555 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl; 2900 2556 } else { 2901 2557 // remove all candidates from the list and then the list itself 2902 CandidateLine.pointlist.clear(); 2903 CandidateLine.pointlist.push_back(Candidate); 2904 Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " 2905 << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl; 2558 class CandidateForTesselation *remover = NULL; 2559 for (CandidateList::iterator it = candidates->begin(); it != candidates->end(); ++it) { 2560 remover = *it; 2561 delete(remover); 2562 } 2563 candidates->clear(); 2564 candidates->push_back(optCandidate); 2565 Log() << Verbose(2) << "ACCEPT: We have found a better candidate: " << *(optCandidate->point) << " with " 2566 << alpha << " and circumsphere's center at " << optCandidate->OptCenter << "." << endl; 2906 2567 } 2907 CandidateLine.ShortestAngle = alpha;2908 Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl;2568 *ShortestAngle = alpha; 2569 //Log() << Verbose(2) << "INFO: There are " << candidates->size() << " candidates in the list now." << endl; 2909 2570 } else { 2910 if (( Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) {2911 Log() << Verbose(1) << "REJECT: Old candidate " << *(Candidate) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl;2571 if ((optCandidate != NULL) && (optCandidate->point != NULL)) { 2572 //Log() << Verbose(2) << "REJECT: Old candidate " << *(optCandidate->point) << " with " << *ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl; 2912 2573 } else { 2913 Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl;2574 //Log() << Verbose(2) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 2914 2575 } 2915 2576 } 2577 2916 2578 } else { 2917 Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl;2579 //Log() << Verbose(2) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl; 2918 2580 } 2919 2581 } else { 2920 Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl;2582 //Log() << Verbose(2) << "REJECT: Three points from " << *BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 2921 2583 } 2922 2584 } else { 2923 2585 if (ThirdNode != NULL) { 2924 Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl;2586 //Log() << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl; 2925 2587 } else { 2926 Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl;2588 //Log() << Verbose(2) << "REJECT: Base triangle " << *BaseLine << " contains Candidate " << *Candidate << "." << endl; 2927 2589 } 2928 2590 } … … 2935 2597 } else { 2936 2598 if (ThirdNode != NULL) 2937 Log() << Verbose( 1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdNode << " is too big!" << endl;2599 Log() << Verbose(2) << "Circumcircle for base line " << *BaseLine << " and third node " << *ThirdNode << " is too big!" << endl; 2938 2600 else 2939 Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl; 2940 } 2941 2942 Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl; 2943 if (CandidateLine.pointlist.size() > 1) { 2944 CandidateLine.pointlist.unique(); 2945 CandidateLine.pointlist.sort(); //SortCandidates); 2946 } 2601 Log() << Verbose(2) << "Circumcircle for base line " << *BaseLine << " is too big!" << endl; 2602 } 2603 2604 //Log() << Verbose(2) << "INFO: Sorting candidate list ..." << endl; 2605 if (candidates->size() > 1) { 2606 candidates->unique(); 2607 candidates->sort(SortCandidates); 2608 } 2609 2610 Log() << Verbose(1) << "End of FindThirdPointForTesselation" << endl; 2947 2611 }; 2948 2612 … … 2954 2618 class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const 2955 2619 { 2956 Info FunctionInfo(__func__);2957 2620 const BoundaryLineSet * lines[2] = { line1, line2 }; 2958 2621 class BoundaryPointSet *node = NULL; … … 2968 2631 { // if insertion fails, we have common endpoint 2969 2632 node = OrderTest.first->second; 2970 Log() << Verbose( 1) << "Common endpoint of lines " << *line12633 Log() << Verbose(5) << "Common endpoint of lines " << *line1 2971 2634 << " and " << *line2 << " is: " << *node << "." << endl; 2972 2635 j = 2; … … 2985 2648 list<BoundaryTriangleSet*> * Tesselation::FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const 2986 2649 { 2987 Info FunctionInfo(__func__);2988 2650 TesselPoint *trianglePoints[3]; 2989 2651 TesselPoint *SecondPoint = NULL; … … 2991 2653 2992 2654 if (LinesOnBoundary.empty()) { 2993 eLog() << Verbose(1) << "Error: There is no tesselation structure to compare the point with, please create one first.";2655 Log() << Verbose(0) << "Error: There is no tesselation structure to compare the point with, please create one first."; 2994 2656 return NULL; 2995 2657 } … … 2999 2661 // check whether closest point is "too close" :), then it's inside 3000 2662 if (trianglePoints[0] == NULL) { 3001 Log() << Verbose( 0) << "Is the only point, no one else is closeby." << endl;2663 Log() << Verbose(2) << "Is the only point, no one else is closeby." << endl; 3002 2664 return NULL; 3003 2665 } 3004 2666 if (trianglePoints[0]->node->DistanceSquared(x) < MYEPSILON) { 3005 Log() << Verbose( 1) << "Point is right on a tesselation point, no nearest triangle." << endl;2667 Log() << Verbose(3) << "Point is right on a tesselation point, no nearest triangle." << endl; 3006 2668 PointMap::const_iterator PointRunner = PointsOnBoundary.find(trianglePoints[0]->nr); 3007 2669 triangles = new list<BoundaryTriangleSet*>; … … 3027 2689 } 3028 2690 } else { 3029 set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(trianglePoints[0]); 3030 TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(connectedPoints, trianglePoints[0], x); 3031 delete(connectedPoints); 2691 list<TesselPoint*> *connectedClosestPoints = GetCircleOfConnectedPoints(trianglePoints[0], x); 3032 2692 if (connectedClosestPoints != NULL) { 3033 2693 trianglePoints[1] = connectedClosestPoints->front(); … … 3037 2697 eLog() << Verbose(1) << "IsInnerPoint encounters serious error, point " << i << " not found." << endl; 3038 2698 } 3039 //Log() << Verbose( 1) << "List of triangle points:" << endl;3040 //Log() << Verbose( 2) << *trianglePoints[i] << endl;2699 //Log() << Verbose(2) << "List of triangle points:" << endl; 2700 //Log() << Verbose(3) << *trianglePoints[i] << endl; 3041 2701 } 3042 2702 3043 2703 triangles = FindTriangles(trianglePoints); 3044 Log() << Verbose( 1) << "List of possible triangles:" << endl;2704 Log() << Verbose(2) << "List of possible triangles:" << endl; 3045 2705 for(list<BoundaryTriangleSet*>::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) 3046 Log() << Verbose( 2) << **Runner << endl;2706 Log() << Verbose(3) << **Runner << endl; 3047 2707 3048 2708 delete(connectedClosestPoints); 3049 2709 } else { 3050 2710 triangles = NULL; 3051 eLog() << Verbose(2) << "There is no circle of connected points!" << endl;2711 Log() << Verbose(1) << "There is no circle of connected points!" << endl; 3052 2712 } 3053 2713 } … … 3069 2729 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const 3070 2730 { 3071 Info FunctionInfo(__func__);3072 2731 class BoundaryTriangleSet *result = NULL; 3073 2732 list<BoundaryTriangleSet*> *triangles = FindClosestTrianglesToPoint(x, LC); … … 3079 2738 if (triangles->size() == 1) { // there is no degenerate case 3080 2739 result = triangles->front(); 3081 Log() << Verbose( 1) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl;2740 Log() << Verbose(2) << "Normal Vector of this triangle is " << result->NormalVector << "." << endl; 3082 2741 } else { 3083 2742 result = triangles->front(); 3084 2743 result->GetCenter(&Center); 3085 2744 Center.SubtractVector(x); 3086 Log() << Verbose( 1) << "Normal Vector of this front side is " << result->NormalVector << "." << endl;2745 Log() << Verbose(2) << "Normal Vector of this front side is " << result->NormalVector << "." << endl; 3087 2746 if (Center.ScalarProduct(&result->NormalVector) < 0) { 3088 2747 result = triangles->back(); 3089 Log() << Verbose( 1) << "Normal Vector of this back side is " << result->NormalVector << "." << endl;2748 Log() << Verbose(2) << "Normal Vector of this back side is " << result->NormalVector << "." << endl; 3090 2749 if (Center.ScalarProduct(&result->NormalVector) < 0) { 3091 2750 eLog() << Verbose(1) << "Front and back side yield NormalVector in wrong direction!" << endl; … … 3106 2765 bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const 3107 2766 { 3108 Info FunctionInfo(__func__);3109 2767 class BoundaryTriangleSet *result = FindClosestTriangleToPoint(&Point, LC); 3110 2768 Vector Center; … … 3116 2774 3117 2775 result->GetCenter(&Center); 3118 Log() << Verbose( 2) << "INFO: Central point of the triangle is " << Center << "." << endl;2776 Log() << Verbose(3) << "INFO: Central point of the triangle is " << Center << "." << endl; 3119 2777 Center.SubtractVector(&Point); 3120 Log() << Verbose( 2) << "INFO: Vector from center to point to test is " << Center << "." << endl;2778 Log() << Verbose(3) << "INFO: Vector from center to point to test is " << Center << "." << endl; 3121 2779 if (Center.ScalarProduct(&result->NormalVector) > -MYEPSILON) { 3122 2780 Log() << Verbose(1) << Point << " is an inner point." << endl; … … 3137 2795 bool Tesselation::IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC) const 3138 2796 { 3139 Info FunctionInfo(__func__);3140 2797 return IsInnerPoint(*(Point->node), LC); 3141 2798 } … … 3149 2806 set<TesselPoint*> * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const 3150 2807 { 3151 Info FunctionInfo(__func__);3152 2808 set<TesselPoint*> *connectedPoints = new set<TesselPoint*>; 3153 2809 class BoundaryPointSet *ReferencePoint = NULL; 3154 2810 TesselPoint* current; 3155 2811 bool takePoint = false; 2812 2813 Log() << Verbose(3) << "Begin of GetAllConnectedPoints" << endl; 3156 2814 3157 2815 // find the respective boundary point … … 3160 2818 ReferencePoint = PointRunner->second; 3161 2819 } else { 3162 eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl;2820 Log() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl; 3163 2821 ReferencePoint = NULL; 3164 2822 } … … 3184 2842 3185 2843 if (takePoint) { 3186 Log() << Verbose( 1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl;2844 Log() << Verbose(5) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl; 3187 2845 connectedPoints->insert(current); 3188 2846 } … … 3196 2854 } 3197 2855 2856 Log() << Verbose(3) << "End of GetAllConnectedPoints" << endl; 3198 2857 return connectedPoints; 3199 2858 }; … … 3207 2866 * 3208 2867 * @param *out output stream for debugging 3209 * @param *SetOfNeighbours all points for which the angle should be calculated3210 2868 * @param *Point of which get all connected points 3211 2869 * @param *Reference Reference vector for zero angle or NULL for no preference 3212 2870 * @return list of the all points linked to the provided one 3213 2871 */ 3214 list<TesselPoint*> * Tesselation::GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const 3215 { 3216 Info FunctionInfo(__func__); 2872 list<TesselPoint*> * Tesselation::GetCircleOfConnectedPoints(const TesselPoint* const Point, const Vector * const Reference) const 2873 { 3217 2874 map<double, TesselPoint*> anglesOfPoints; 2875 set<TesselPoint*> *connectedPoints = GetAllConnectedPoints(Point); 3218 2876 list<TesselPoint*> *connectedCircle = new list<TesselPoint*>; 3219 2877 Vector center; … … 3223 2881 Vector helper; 3224 2882 3225 if ( SetOfNeighbours == NULL) {3226 eLog() << Verbose(2) << "Could not find any connected points!" << endl;2883 if (connectedPoints == NULL) { 2884 Log() << Verbose(2) << "Could not find any connected points!" << endl; 3227 2885 delete(connectedCircle); 3228 2886 return NULL; 3229 2887 } 2888 Log() << Verbose(2) << "Begin of GetCircleOfConnectedPoints" << endl; 3230 2889 3231 2890 // calculate central point 3232 for (set<TesselPoint*>::const_iterator TesselRunner = SetOfNeighbours->begin(); TesselRunner != SetOfNeighbours->end(); TesselRunner++)2891 for (set<TesselPoint*>::const_iterator TesselRunner = connectedPoints->begin(); TesselRunner != connectedPoints->end(); TesselRunner++) 3233 2892 center.AddVector((*TesselRunner)->node); 3234 2893 //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size() 3235 2894 // << "; scale factor " << 1.0/connectedPoints.size(); 3236 center.Scale(1.0/ SetOfNeighbours->size());3237 Log() << Verbose( 1) << "INFO: Calculated center of all circle points is " << center << "." << endl;2895 center.Scale(1.0/connectedPoints->size()); 2896 Log() << Verbose(4) << "INFO: Calculated center of all circle points is " << center << "." << endl; 3238 2897 3239 2898 // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points … … 3241 2900 PlaneNormal.SubtractVector(¢er); 3242 2901 PlaneNormal.Normalize(); 3243 Log() << Verbose( 1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl;2902 Log() << Verbose(4) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl; 3244 2903 3245 2904 // construct one orthogonal vector … … 3250 2909 } 3251 2910 if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) { 3252 Log() << Verbose( 1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl;3253 AngleZero.CopyVector((* SetOfNeighbours->begin())->node);2911 Log() << Verbose(4) << "Using alternatively " << *(*connectedPoints->begin())->node << " as angle 0 referencer." << endl; 2912 AngleZero.CopyVector((*connectedPoints->begin())->node); 3254 2913 AngleZero.SubtractVector(Point->node); 3255 2914 AngleZero.ProjectOntoPlane(&PlaneNormal); … … 3259 2918 } 3260 2919 } 3261 Log() << Verbose( 1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl;2920 Log() << Verbose(4) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl; 3262 2921 if (AngleZero.NormSquared() > MYEPSILON) 3263 2922 OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero); 3264 2923 else 3265 2924 OrthogonalVector.MakeNormalVector(&PlaneNormal); 3266 Log() << Verbose( 1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl;2925 Log() << Verbose(4) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl; 3267 2926 3268 2927 // go through all connected points and calculate angle 3269 for (set<TesselPoint*>::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) {2928 for (set<TesselPoint*>::iterator listRunner = connectedPoints->begin(); listRunner != connectedPoints->end(); listRunner++) { 3270 2929 helper.CopyVector((*listRunner)->node); 3271 2930 helper.SubtractVector(Point->node); 3272 2931 helper.ProjectOntoPlane(&PlaneNormal); 3273 2932 double angle = GetAngle(helper, AngleZero, OrthogonalVector); 3274 Log() << Verbose( 0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl;2933 Log() << Verbose(3) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl; 3275 2934 anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner))); 3276 2935 } … … 3279 2938 connectedCircle->push_back(AngleRunner->second); 3280 2939 } 2940 2941 delete(connectedPoints); 2942 2943 Log() << Verbose(2) << "End of GetCircleOfConnectedPoints" << endl; 3281 2944 3282 2945 return connectedCircle; … … 3291 2954 list<list<TesselPoint*> *> * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const 3292 2955 { 3293 Info FunctionInfo(__func__);3294 2956 map<double, TesselPoint*> anglesOfPoints; 3295 2957 list<list<TesselPoint*> *> *ListOfPaths = new list<list<TesselPoint*> *>; … … 3336 2998 StartLine = CurrentLine; 3337 2999 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint); 3338 Log() << Verbose( 1)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl;3000 Log() << Verbose(3)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl; 3339 3001 do { 3340 3002 // push current one 3341 Log() << Verbose( 1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;3003 Log() << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl; 3342 3004 connectedPath->push_back(CurrentPoint->node); 3343 3005 3344 3006 // find next triangle 3345 3007 for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) { 3346 Log() << Verbose( 1) << "INFO: Inspecting triangle " << *Runner->second << "." << endl;3008 Log() << Verbose(3) << "INFO: Inspecting triangle " << *Runner->second << "." << endl; 3347 3009 if ((Runner->second != triangle)) { // look for first triangle not equal to old one 3348 3010 triangle = Runner->second; … … 3351 3013 if (!TriangleRunner->second) { 3352 3014 TriangleRunner->second = true; 3353 Log() << Verbose( 1) << "INFO: Connecting triangle is " << *triangle << "." << endl;3015 Log() << Verbose(3) << "INFO: Connecting triangle is " << *triangle << "." << endl; 3354 3016 break; 3355 3017 } else { 3356 Log() << Verbose( 1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl;3018 Log() << Verbose(3) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl; 3357 3019 triangle = NULL; 3358 3020 } … … 3369 3031 if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point 3370 3032 CurrentLine = triangle->lines[i]; 3371 Log() << Verbose( 1) << "INFO: Connecting line is " << *CurrentLine << "." << endl;3033 Log() << Verbose(3) << "INFO: Connecting line is " << *CurrentLine << "." << endl; 3372 3034 break; 3373 3035 } … … 3383 3045 } while (CurrentLine != StartLine); 3384 3046 // last point is missing, as it's on start line 3385 Log() << Verbose( 1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl;3047 Log() << Verbose(3) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl; 3386 3048 if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back()) 3387 3049 connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node); … … 3389 3051 ListOfPaths->push_back(connectedPath); 3390 3052 } else { 3391 Log() << Verbose( 1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl;3053 Log() << Verbose(3) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl; 3392 3054 } 3393 3055 } … … 3407 3069 list<list<TesselPoint*> *> * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const 3408 3070 { 3409 Info FunctionInfo(__func__);3410 3071 list<list<TesselPoint*> *> *ListofPaths = GetPathsOfConnectedPoints(Point); 3411 3072 list<list<TesselPoint*> *> *ListofClosedPaths = new list<list<TesselPoint*> *>; … … 3421 3082 connectedPath = *ListRunner; 3422 3083 3423 Log() << Verbose( 1) << "INFO: Current path is " << connectedPath << "." << endl;3084 Log() << Verbose(2) << "INFO: Current path is " << connectedPath << "." << endl; 3424 3085 3425 3086 // go through list, look for reappearance of starting Point and count … … 3431 3092 if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point 3432 3093 // we have a closed circle from Marker to new Marker 3433 Log() << Verbose( 1) << count+1 << ". closed path consists of: ";3094 Log() << Verbose(3) << count+1 << ". closed path consists of: "; 3434 3095 newPath = new list<TesselPoint*>; 3435 3096 list<TesselPoint*>::iterator CircleSprinter = Marker; … … 3447 3108 } 3448 3109 } 3449 Log() << Verbose( 1) << "INFO: " << count << " closed additional path(s) have been created." << endl;3110 Log() << Verbose(3) << "INFO: " << count << " closed additional path(s) have been created." << endl; 3450 3111 3451 3112 // delete list of paths … … 3469 3130 set<BoundaryTriangleSet*> *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const 3470 3131 { 3471 Info FunctionInfo(__func__);3472 3132 set<BoundaryTriangleSet*> *connectedTriangles = new set<BoundaryTriangleSet*>; 3473 3133 … … 3508 3168 return 0.; 3509 3169 } else 3510 Log() << Verbose( 0) << "Removing point " << *point << " from tesselated boundary ..." << endl;3170 Log() << Verbose(2) << "Removing point " << *point << " from tesselated boundary ..." << endl; 3511 3171 3512 3172 // copy old location for the volume … … 3538 3198 NormalVector.Zero(); 3539 3199 for (map<class BoundaryTriangleSet *, int>::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) { 3540 Log() << Verbose( 1) << "INFO: Removing triangle " << *(Runner->first) << "." << endl;3200 Log() << Verbose(3) << "INFO: Removing triangle " << *(Runner->first) << "." << endl; 3541 3201 NormalVector.SubtractVector(&Runner->first->NormalVector); // has to point inward 3542 3202 RemoveTesselationTriangle(Runner->first); … … 3568 3228 smallestangle = 0.; 3569 3229 for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) { 3570 Log() << Verbose( 1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;3230 Log() << Verbose(3) << "INFO: MiddleNode is " << **MiddleNode << "." << endl; 3571 3231 // construct vectors to next and previous neighbour 3572 3232 StartNode = MiddleNode; … … 3596 3256 MiddleNode = EndNode; 3597 3257 if (MiddleNode == connectedPath->end()) { 3598 eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl;3599 performCriticalExit();3258 Log() << Verbose(1) << "CRITICAL: Could not find a smallest angle!" << endl; 3259 exit(255); 3600 3260 } 3601 3261 StartNode = MiddleNode; … … 3606 3266 if (EndNode == connectedPath->end()) 3607 3267 EndNode = connectedPath->begin(); 3608 Log() << Verbose( 2) << "INFO: StartNode is " << **StartNode << "." << endl;3609 Log() << Verbose( 2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl;3610 Log() << Verbose( 2) << "INFO: EndNode is " << **EndNode << "." << endl;3611 Log() << Verbose( 1) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl;3268 Log() << Verbose(4) << "INFO: StartNode is " << **StartNode << "." << endl; 3269 Log() << Verbose(4) << "INFO: MiddleNode is " << **MiddleNode << "." << endl; 3270 Log() << Verbose(4) << "INFO: EndNode is " << **EndNode << "." << endl; 3271 Log() << Verbose(3) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl; 3612 3272 TriangleCandidates[0] = *StartNode; 3613 3273 TriangleCandidates[1] = *MiddleNode; … … 3615 3275 triangle = GetPresentTriangle(TriangleCandidates); 3616 3276 if (triangle != NULL) { 3617 eLog() << Verbose( 0) << "New triangle already present, skipping!" << endl;3277 eLog() << Verbose(2) << "New triangle already present, skipping!" << endl; 3618 3278 StartNode++; 3619 3279 MiddleNode++; … … 3627 3287 continue; 3628 3288 } 3629 Log() << Verbose( 3) << "Adding new triangle points."<< endl;3289 Log() << Verbose(5) << "Adding new triangle points."<< endl; 3630 3290 AddTesselationPoint(*StartNode, 0); 3631 3291 AddTesselationPoint(*MiddleNode, 1); 3632 3292 AddTesselationPoint(*EndNode, 2); 3633 Log() << Verbose( 3) << "Adding new triangle lines."<< endl;3293 Log() << Verbose(5) << "Adding new triangle lines."<< endl; 3634 3294 AddTesselationLine(TPS[0], TPS[1], 0); 3635 3295 AddTesselationLine(TPS[0], TPS[2], 1); … … 3646 3306 // prepare nodes for next triangle 3647 3307 StartNode = EndNode; 3648 Log() << Verbose( 2) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl;3308 Log() << Verbose(4) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl; 3649 3309 connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles) 3650 3310 if (connectedPath->size() == 2) { // we are done … … 3653 3313 break; 3654 3314 } else if (connectedPath->size() < 2) { // something's gone wrong! 3655 eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl;3656 performCriticalExit();3315 Log() << Verbose(1) << "CRITICAL: There are only two endpoints left!" << endl; 3316 exit(255); 3657 3317 } else { 3658 3318 MiddleNode = StartNode; … … 3682 3342 if (maxgain != 0) { 3683 3343 volume += maxgain; 3684 Log() << Verbose( 1) << "Flipping baseline with highest volume" << **Candidate << "." << endl;3344 Log() << Verbose(3) << "Flipping baseline with highest volume" << **Candidate << "." << endl; 3685 3345 OtherBase = FlipBaseline(*Candidate); 3686 3346 NewLines.erase(Candidate); … … 3693 3353 delete(connectedPath); 3694 3354 } 3695 Log() << Verbose( 0) << count << " triangles were created." << endl;3355 Log() << Verbose(1) << count << " triangles were created." << endl; 3696 3356 } else { 3697 3357 while (!ListOfClosedPaths->empty()) { … … 3701 3361 delete(connectedPath); 3702 3362 } 3703 Log() << Verbose( 0) << "No need to create any triangles." << endl;3363 Log() << Verbose(1) << "No need to create any triangles." << endl; 3704 3364 } 3705 3365 delete(ListOfClosedPaths); 3706 3366 3707 Log() << Verbose( 0) << "Removed volume is " << volume << "." << endl;3367 Log() << Verbose(1) << "Removed volume is " << volume << "." << endl; 3708 3368 3709 3369 return volume; … … 3722 3382 list<BoundaryTriangleSet*> *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const 3723 3383 { 3724 Info FunctionInfo(__func__);3725 3384 list<BoundaryTriangleSet*> *result = new list<BoundaryTriangleSet*>; 3726 3385 LineMap::const_iterator FindLine; … … 3763 3422 } 3764 3423 3765 struct BoundaryLineSetCompare {3766 bool operator() (const BoundaryLineSet * const a, const BoundaryLineSet * const b) {3767 int lowerNra = -1;3768 int lowerNrb = -1;3769 3770 if (a->endpoints[0] < a->endpoints[1])3771 lowerNra = 0;3772 else3773 lowerNra = 1;3774 3775 if (b->endpoints[0] < b->endpoints[1])3776 lowerNrb = 0;3777 else3778 lowerNrb = 1;3779 3780 if (a->endpoints[lowerNra] < b->endpoints[lowerNrb])3781 return true;3782 else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb])3783 return false;3784 else { // both lower-numbered endpoints are the same ...3785 if (a->endpoints[(lowerNra+1)%2] < b->endpoints[(lowerNrb+1)%2])3786 return true;3787 else if (a->endpoints[(lowerNra+1)%2] > b->endpoints[(lowerNrb+1)%2])3788 return false;3789 }3790 return false;3791 };3792 };3793 3794 #define UniqueLines set < class BoundaryLineSet *, BoundaryLineSetCompare>3795 3796 3424 /** 3797 3425 * Finds all degenerated lines within the tesselation structure. … … 3802 3430 map<int, int> * Tesselation::FindAllDegeneratedLines() 3803 3431 { 3804 Info FunctionInfo(__func__); 3805 UniqueLines AllLines; 3432 map<int, class BoundaryLineSet *> AllLines; 3806 3433 map<int, int> * DegeneratedLines = new map<int, int>; 3807 3434 3808 3435 // sanity check 3809 3436 if (LinesOnBoundary.empty()) { 3810 eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.";3437 Log() << Verbose(1) << "Warning: FindAllDegeneratedTriangles() was called without any tesselation structure."; 3811 3438 return DegeneratedLines; 3812 3439 } 3813 3440 3814 3441 LineMap::iterator LineRunner1; 3815 pair< UniqueLines::iterator, bool> tester;3442 pair<LineMap::iterator, bool> tester; 3816 3443 for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) { 3817 tester = AllLines.insert( LineRunner1->second);3818 if ( !tester.second) { // found degenerated line3819 DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr) );3820 DegeneratedLines->insert ( pair<int, int> ( (*tester.first)->Nr, LineRunner1->second->Nr) );3444 tester = AllLines.insert( pair<int,BoundaryLineSet *> (LineRunner1->second->endpoints[0]->Nr, LineRunner1->second) ); 3445 if ((!tester.second) && (tester.first->second->endpoints[1]->Nr == LineRunner1->second->endpoints[1]->Nr)) { // found degenerated line 3446 DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, tester.first->second->Nr) ); 3447 DegeneratedLines->insert ( pair<int, int> (tester.first->second->Nr, LineRunner1->second->Nr) ); 3821 3448 } 3822 3449 } … … 3824 3451 AllLines.clear(); 3825 3452 3826 Log() << Verbose( 0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl;3453 Log() << Verbose(1) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl; 3827 3454 map<int,int>::iterator it; 3828 for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) { 3829 const LineMap::const_iterator Line1 = LinesOnBoundary.find((*it).first); 3830 const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second); 3831 if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end()) 3832 Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl; 3833 else 3834 eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl; 3835 } 3455 for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) 3456 Log() << Verbose(2) << (*it).first << " => " << (*it).second << endl; 3836 3457 3837 3458 return DegeneratedLines; … … 3846 3467 map<int, int> * Tesselation::FindAllDegeneratedTriangles() 3847 3468 { 3848 Info FunctionInfo(__func__);3849 3469 map<int, int> * DegeneratedLines = FindAllDegeneratedLines(); 3850 3470 map<int, int> * DegeneratedTriangles = new map<int, int>; … … 3874 3494 delete(DegeneratedLines); 3875 3495 3876 Log() << Verbose( 0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl;3496 Log() << Verbose(1) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl; 3877 3497 map<int,int>::iterator it; 3878 3498 for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++) 3879 Log() << Verbose( 0) << (*it).first << " => " << (*it).second << endl;3499 Log() << Verbose(2) << (*it).first << " => " << (*it).second << endl; 3880 3500 3881 3501 return DegeneratedTriangles; … … 3888 3508 void Tesselation::RemoveDegeneratedTriangles() 3889 3509 { 3890 Info FunctionInfo(__func__);3891 3510 map<int, int> * DegeneratedTriangles = FindAllDegeneratedTriangles(); 3892 3511 TriangleMap::iterator finder; 3893 3512 BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL; 3894 3513 int count = 0; 3514 3515 Log() << Verbose(1) << "Begin of RemoveDegeneratedTriangles" << endl; 3895 3516 3896 3517 for (map<int, int>::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); … … 3951 3572 // erase the pair 3952 3573 count += (int) DegeneratedTriangles->erase(triangle->Nr); 3953 Log() << Verbose( 0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl;3574 Log() << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl; 3954 3575 RemoveTesselationTriangle(triangle); 3955 3576 count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr); 3956 Log() << Verbose( 0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl;3577 Log() << Verbose(1) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl; 3957 3578 RemoveTesselationTriangle(partnerTriangle); 3958 3579 } else { 3959 Log() << Verbose( 0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle3580 Log() << Verbose(1) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle 3960 3581 << " and its partner " << *partnerTriangle << " because it is essential for at" 3961 3582 << " least one of the endpoints to be kept in the tesselation structure." << endl; … … 3963 3584 } 3964 3585 delete(DegeneratedTriangles); 3965 if (count > 0) 3966 LastTriangle = NULL; 3967 3968 Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl; 3586 3587 Log() << Verbose(1) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl; 3588 Log() << Verbose(1) << "End of RemoveDegeneratedTriangles" << endl; 3969 3589 } 3970 3590 … … 3979 3599 void Tesselation::AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC) 3980 3600 { 3981 Info FunctionInfo(__func__); 3601 Log() << Verbose(2) << "Begin of AddBoundaryPointByDegeneratedTriangle" << endl; 3602 3982 3603 // find nearest boundary point 3983 3604 class TesselPoint *BackupPoint = NULL; … … 3995 3616 return; 3996 3617 } 3997 Log() << Verbose( 0) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl;3618 Log() << Verbose(2) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl; 3998 3619 3999 3620 // go through its lines and find the best one to split … … 4028 3649 4029 3650 // create new triangle to connect point (connects automatically with the missing spot of the chosen line) 4030 Log() << Verbose( 2) << "Adding new triangle points."<< endl;3651 Log() << Verbose(5) << "Adding new triangle points."<< endl; 4031 3652 AddTesselationPoint((BestLine->endpoints[0]->node), 0); 4032 3653 AddTesselationPoint((BestLine->endpoints[1]->node), 1); 4033 3654 AddTesselationPoint(point, 2); 4034 Log() << Verbose( 2) << "Adding new triangle lines."<< endl;3655 Log() << Verbose(5) << "Adding new triangle lines."<< endl; 4035 3656 AddTesselationLine(TPS[0], TPS[1], 0); 4036 3657 AddTesselationLine(TPS[0], TPS[2], 1); … … 4039 3660 BTS->GetNormalVector(TempTriangle->NormalVector); 4040 3661 BTS->NormalVector.Scale(-1.); 4041 Log() << Verbose( 1) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl;3662 Log() << Verbose(3) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl; 4042 3663 AddTesselationTriangle(); 4043 3664 4044 3665 // create other side of this triangle and close both new sides of the first created triangle 4045 Log() << Verbose( 2) << "Adding new triangle points."<< endl;3666 Log() << Verbose(5) << "Adding new triangle points."<< endl; 4046 3667 AddTesselationPoint((BestLine->endpoints[0]->node), 0); 4047 3668 AddTesselationPoint((BestLine->endpoints[1]->node), 1); 4048 3669 AddTesselationPoint(point, 2); 4049 Log() << Verbose( 2) << "Adding new triangle lines."<< endl;3670 Log() << Verbose(5) << "Adding new triangle lines."<< endl; 4050 3671 AddTesselationLine(TPS[0], TPS[1], 0); 4051 3672 AddTesselationLine(TPS[0], TPS[2], 1); … … 4053 3674 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 4054 3675 BTS->GetNormalVector(TempTriangle->NormalVector); 4055 Log() << Verbose( 1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl;3676 Log() << Verbose(3) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl; 4056 3677 AddTesselationTriangle(); 4057 3678 … … 4060 3681 if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) { 4061 3682 if (BestLine == BTS->lines[i]){ 4062 eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl;4063 performCriticalExit();3683 Log() << Verbose(1) << "CRITICAL: BestLine is same as found line, something's wrong here!" << endl; 3684 exit(255); 4064 3685 } 4065 3686 BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) ); … … 4068 3689 } 4069 3690 } 3691 3692 // exit 3693 Log() << Verbose(2) << "End of AddBoundaryPointByDegeneratedTriangle" << endl; 4070 3694 }; 4071 3695 … … 4077 3701 void Tesselation::Output(const char *filename, const PointCloud * const cloud) 4078 3702 { 4079 Info FunctionInfo(__func__);4080 3703 ofstream *tempstream = NULL; 4081 3704 string NameofTempFile; … … 4090 3713 NameofTempFile.erase(npos, 1); 4091 3714 NameofTempFile.append(TecplotSuffix); 4092 Log() << Verbose( 0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";3715 Log() << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 4093 3716 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 4094 3717 WriteTecplotFile(tempstream, this, cloud, TriangleFilesWritten); … … 4104 3727 NameofTempFile.erase(npos, 1); 4105 3728 NameofTempFile.append(Raster3DSuffix); 4106 Log() << Verbose( 0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n";3729 Log() << Verbose(1) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 4107 3730 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 4108 3731 WriteRaster3dFile(tempstream, this, cloud); … … 4116 3739 TriangleFilesWritten++; 4117 3740 }; 4118 4119 struct BoundaryPolygonSetCompare {4120 bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const {4121 if (s1->endpoints.size() < s2->endpoints.size())4122 return true;4123 else if (s1->endpoints.size() > s2->endpoints.size())4124 return false;4125 else { // equality of number of endpoints4126 PointSet::const_iterator Walker1 = s1->endpoints.begin();4127 PointSet::const_iterator Walker2 = s2->endpoints.begin();4128 while ((Walker1 != s1->endpoints.end()) || (Walker2 != s2->endpoints.end())) {4129 if ((*Walker1)->Nr < (*Walker2)->Nr)4130 return true;4131 else if ((*Walker1)->Nr > (*Walker2)->Nr)4132 return false;4133 Walker1++;4134 Walker2++;4135 }4136 return false;4137 }4138 }4139 };4140 4141 #define UniquePolygonSet set < BoundaryPolygonSet *, BoundaryPolygonSetCompare>4142 4143 /** Finds all degenerated polygons and calls ReTesselateDegeneratedPolygon()/4144 * \return number of polygons found4145 */4146 int Tesselation::CorrectAllDegeneratedPolygons()4147 {4148 Info FunctionInfo(__func__);4149 4150 /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector4151 map <int, int> *DegeneratedTriangles = FindAllDegeneratedTriangles();4152 set < BoundaryPointSet *> EndpointCandidateList;4153 pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester;4154 pair < map < int, Vector *>::iterator, bool > TriangleInsertionTester;4155 for (PointMap::const_iterator Runner = PointsOnBoundary.begin(); Runner != PointsOnBoundary.end(); Runner++) {4156 Log() << Verbose(0) << "Current point is " << *Runner->second << "." << endl;4157 map < int, Vector *> TriangleVectors;4158 // gather all NormalVectors4159 Log() << Verbose(1) << "Gathering triangles ..." << endl;4160 for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++)4161 for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) {4162 if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) {4163 TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) );4164 if (TriangleInsertionTester.second)4165 Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl;4166 } else {4167 Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl;4168 }4169 }4170 // check whether there are two that are parallel4171 Log() << Verbose(1) << "Finding two parallel triangles ..." << endl;4172 for (map < int, Vector *>::iterator VectorWalker = TriangleVectors.begin(); VectorWalker != TriangleVectors.end(); VectorWalker++)4173 for (map < int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++)4174 if (VectorWalker != VectorRunner) { // skip equals4175 const double SCP = VectorWalker->second->ScalarProduct(VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles4176 Log() << Verbose(1) << "Checking " << *VectorWalker->second<< " against " << *VectorRunner->second << ": " << SCP << endl;4177 if (fabs(SCP + 1.) < ParallelEpsilon) {4178 InsertionTester = EndpointCandidateList.insert((Runner->second));4179 if (InsertionTester.second)4180 Log() << Verbose(0) << " Adding " << *Runner->second << " to endpoint candidate list." << endl;4181 // and break out of both loops4182 VectorWalker = TriangleVectors.end();4183 VectorRunner = TriangleVectors.end();4184 break;4185 }4186 }4187 }4188 4189 /// 3. Find connected endpoint candidates and put them into a polygon4190 UniquePolygonSet ListofDegeneratedPolygons;4191 BoundaryPointSet *Walker = NULL;4192 BoundaryPointSet *OtherWalker = NULL;4193 BoundaryPolygonSet *Current = NULL;4194 stack <BoundaryPointSet*> ToCheckConnecteds;4195 while (!EndpointCandidateList.empty()) {4196 Walker = *(EndpointCandidateList.begin());4197 if (Current == NULL) { // create a new polygon with current candidate4198 Log() << Verbose(0) << "Starting new polygon set at point " << *Walker << endl;4199 Current = new BoundaryPolygonSet;4200 Current->endpoints.insert(Walker);4201 EndpointCandidateList.erase(Walker);4202 ToCheckConnecteds.push(Walker);4203 }4204 4205 // go through to-check stack4206 while (!ToCheckConnecteds.empty()) {4207 Walker = ToCheckConnecteds.top(); // fetch ...4208 ToCheckConnecteds.pop(); // ... and remove4209 for (LineMap::const_iterator LineWalker = Walker->lines.begin(); LineWalker != Walker->lines.end(); LineWalker++) {4210 OtherWalker = (LineWalker->second)->GetOtherEndpoint(Walker);4211 Log() << Verbose(1) << "Checking " << *OtherWalker << endl;4212 set < BoundaryPointSet *>::iterator Finder = EndpointCandidateList.find(OtherWalker);4213 if (Finder != EndpointCandidateList.end()) { // found a connected partner4214 Log() << Verbose(1) << " Adding to polygon." << endl;4215 Current->endpoints.insert(OtherWalker);4216 EndpointCandidateList.erase(Finder); // remove from candidates4217 ToCheckConnecteds.push(OtherWalker); // but check its partners too4218 } else {4219 Log() << Verbose(1) << " is not connected to " << *Walker << endl;4220 }4221 }4222 }4223 4224 Log() << Verbose(0) << "Final polygon is " << *Current << endl;4225 ListofDegeneratedPolygons.insert(Current);4226 Current = NULL;4227 }4228 4229 const int counter = ListofDegeneratedPolygons.size();4230 4231 Log() << Verbose(0) << "The following " << counter << " degenerated polygons have been found: " << endl;4232 for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++)4233 Log() << Verbose(0) << " " << **PolygonRunner << endl;4234 4235 /// 4. Go through all these degenerated polygons4236 for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++) {4237 stack <int> TriangleNrs;4238 Vector NormalVector;4239 /// 4a. Gather all triangles of this polygon4240 TriangleSet *T = (*PolygonRunner)->GetAllContainedTrianglesFromEndpoints();4241 4242 // check whether number is bigger than 2, otherwise it's just a simply degenerated one and nothing to do.4243 if (T->size() == 2) {4244 Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl;4245 delete(T);4246 continue;4247 }4248 4249 // check whether number is even4250 // If this case occurs, we have to think about it!4251 // The Problem is probably due to two degenerated polygons being connected by a bridging, non-degenerated polygon, as somehow one node has4252 // connections to either polygon ...4253 if (T->size() % 2 != 0) {4254 eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl;4255 performCriticalExit();4256 }4257 4258 TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator4259 /// 4a. Get NormalVector for one side (this is "front")4260 NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);4261 Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl;4262 TriangleWalker++;4263 TriangleSet::iterator TriangleSprinter = TriangleWalker; // is the inner advanced iterator4264 /// 4b. Remove all triangles whose NormalVector is in opposite direction (i.e. "back")4265 BoundaryTriangleSet *triangle = NULL;4266 while (TriangleSprinter != T->end()) {4267 TriangleWalker = TriangleSprinter;4268 triangle = *TriangleWalker;4269 TriangleSprinter++;4270 Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl;4271 if (triangle->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list4272 Log() << Verbose(1) << " Removing ... " << endl;4273 TriangleNrs.push(triangle->Nr);4274 T->erase(TriangleWalker);4275 RemoveTesselationTriangle(triangle);4276 } else4277 Log() << Verbose(1) << " Keeping ... " << endl;4278 }4279 /// 4c. Copy all "front" triangles but with inverse NormalVector4280 TriangleWalker = T->begin();4281 while (TriangleWalker != T->end()) { // go through all front triangles4282 Log() << Verbose(1) << " Re-creating triangle " << **TriangleWalker << " with NormalVector " << (*TriangleWalker)->NormalVector << endl;4283 for (int i = 0; i < 3; i++)4284 AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i);4285 AddTesselationLine(TPS[0], TPS[1], 0);4286 AddTesselationLine(TPS[0], TPS[2], 1);4287 AddTesselationLine(TPS[1], TPS[2], 2);4288 if (TriangleNrs.empty())4289 eLog() << Verbose(0) << "No more free triangle numbers!" << endl;4290 BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ...4291 AddTesselationTriangle(); // ... and add4292 TriangleNrs.pop();4293 BTS->NormalVector.CopyVector(&(*TriangleWalker)->NormalVector);4294 BTS->NormalVector.Scale(-1.);4295 TriangleWalker++;4296 }4297 if (!TriangleNrs.empty()) {4298 eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl;4299 }4300 delete(T); // remove the triangleset4301 }4302 4303 map<int, int> * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles();4304 Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl;4305 map<int,int>::iterator it;4306 for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++)4307 Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl;4308 delete(SimplyDegeneratedTriangles);4309 4310 /// 5. exit4311 UniquePolygonSet::iterator PolygonRunner;4312 while (!ListofDegeneratedPolygons.empty()) {4313 PolygonRunner = ListofDegeneratedPolygons.begin();4314 delete(*PolygonRunner);4315 ListofDegeneratedPolygons.erase(PolygonRunner);4316 }4317 4318 return counter;4319 };
Note:
See TracChangeset
for help on using the changeset viewer.