Changes in src/tesselation.cpp [b1a6d8:244a84]
- File:
-
- 1 edited
-
src/tesselation.cpp (modified) (257 diffs)
Legend:
- Unmodified
- Added
- Removed
-
src/tesselation.cpp
rb1a6d8 r244a84 7 7 8 8 #include <fstream> 9 #include <assert.h>10 9 11 10 #include "helpers.hpp" … … 15 14 #include "tesselation.hpp" 16 15 #include "tesselationhelpers.hpp" 17 #include "triangleintersectionlist.hpp"18 16 #include "vector.hpp" 19 17 #include "verbose.hpp" … … 26 24 */ 27 25 BoundaryPointSet::BoundaryPointSet() : 28 LinesCount(0), value(0.), Nr(-1) 29 { 30 Info FunctionInfo(__func__); 31 DoLog(1) && (Log() << Verbose(1) << "Adding noname." << endl); 32 } 33 ; 26 LinesCount(0), 27 value(0.), 28 Nr(-1) 29 { 30 Info FunctionInfo(__func__); 31 Log() << Verbose(1) << "Adding noname." << endl; 32 }; 34 33 35 34 /** Constructor of BoundaryPointSet with Tesselpoint. … … 37 36 */ 38 37 BoundaryPointSet::BoundaryPointSet(TesselPoint * const Walker) : 39 LinesCount(0), node(Walker), value(0.), Nr(Walker->nr) 40 { 41 Info FunctionInfo(__func__); 42 DoLog(1) && (Log() << Verbose(1) << "Adding Node " << *Walker << endl); 43 } 44 ; 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; 45 }; 45 46 46 47 /** Destructor of BoundaryPointSet. … … 50 51 BoundaryPointSet::~BoundaryPointSet() 51 52 { 52 Info FunctionInfo(__func__);53 Info FunctionInfo(__func__); 53 54 //Log() << Verbose(0) << "Erasing point nr. " << Nr << "." << endl; 54 55 if (!lines.empty()) 55 DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl);56 eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some lines." << endl; 56 57 node = NULL; 57 } 58 ; 58 }; 59 59 60 60 /** Add a line to the LineMap of this point. … … 63 63 void BoundaryPointSet::AddLine(BoundaryLineSet * const line) 64 64 { 65 Info FunctionInfo(__func__); 66 DoLog(1) && (Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." << endl); 67 if (line->endpoints[0] == this) { 68 lines.insert(LinePair(line->endpoints[1]->Nr, line)); 69 } else { 70 lines.insert(LinePair(line->endpoints[0]->Nr, line)); 71 } 65 Info FunctionInfo(__func__); 66 Log() << Verbose(1) << "Adding " << *this << " to line " << *line << "." 67 << endl; 68 if (line->endpoints[0] == this) 69 { 70 lines.insert(LinePair(line->endpoints[1]->Nr, line)); 71 } 72 else 73 { 74 lines.insert(LinePair(line->endpoints[0]->Nr, line)); 75 } 72 76 LinesCount++; 73 } 74 ; 77 }; 75 78 76 79 /** output operator for BoundaryPointSet. … … 90 93 */ 91 94 BoundaryLineSet::BoundaryLineSet() : 92 Nr(-1)93 { 94 Info FunctionInfo(__func__);95 Nr(-1) 96 { 97 Info FunctionInfo(__func__); 95 98 for (int i = 0; i < 2; i++) 96 99 endpoints[i] = NULL; 97 } 98 ; 100 }; 99 101 100 102 /** Constructor of BoundaryLineSet with two endpoints. … … 105 107 BoundaryLineSet::BoundaryLineSet(BoundaryPointSet * const Point[2], const int number) 106 108 { 107 Info FunctionInfo(__func__);109 Info FunctionInfo(__func__); 108 110 // set number 109 111 Nr = number; … … 116 118 skipped = false; 117 119 // clear triangles list 118 DoLog(0) && (Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl); 119 } 120 ; 120 Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl; 121 }; 121 122 122 123 /** Constructor of BoundaryLineSet with two endpoints. … … 139 140 skipped = false; 140 141 // clear triangles list 141 DoLog(0) && (Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl); 142 } 143 ; 142 Log() << Verbose(0) << "New Line with endpoints " << *this << "." << endl; 143 }; 144 144 145 145 /** Destructor for BoundaryLineSet. … … 149 149 BoundaryLineSet::~BoundaryLineSet() 150 150 { 151 Info FunctionInfo(__func__);151 Info FunctionInfo(__func__); 152 152 int Numbers[2]; 153 153 … … 180 180 //Log() << Verbose(0) << *endpoints[i] << " has no more lines it's attached to, erasing." << endl; 181 181 if (endpoints[i] != NULL) { 182 delete (endpoints[i]);182 delete(endpoints[i]); 183 183 endpoints[i] = NULL; 184 184 } … … 187 187 } 188 188 if (!triangles.empty()) 189 DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl); 190 } 191 ; 189 eLog() << Verbose(2) << "Memory Leak! I " << *this << " am still connected to some triangles." << endl; 190 }; 192 191 193 192 /** Add triangle to TriangleMap of this boundary line. … … 196 195 void BoundaryLineSet::AddTriangle(BoundaryTriangleSet * const triangle) 197 196 { 198 Info FunctionInfo(__func__);199 DoLog(0) && (Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl);197 Info FunctionInfo(__func__); 198 Log() << Verbose(0) << "Add " << triangle->Nr << " to line " << *this << "." << endl; 200 199 triangles.insert(TrianglePair(triangle->Nr, triangle)); 201 } 202 ; 200 }; 203 201 204 202 /** Checks whether we have a common endpoint with given \a *line. … … 208 206 bool BoundaryLineSet::IsConnectedTo(const BoundaryLineSet * const line) const 209 207 { 210 Info FunctionInfo(__func__);208 Info FunctionInfo(__func__); 211 209 if ((endpoints[0] == line->endpoints[0]) || (endpoints[1] == line->endpoints[0]) || (endpoints[0] == line->endpoints[1]) || (endpoints[1] == line->endpoints[1])) 212 210 return true; 213 211 else 214 212 return false; 215 } 216 ; 213 }; 217 214 218 215 /** Checks whether the adjacent triangles of a baseline are convex or not. … … 224 221 bool BoundaryLineSet::CheckConvexityCriterion() const 225 222 { 226 Info FunctionInfo(__func__);223 Info FunctionInfo(__func__); 227 224 Vector BaseLineCenter, BaseLineNormal, BaseLine, helper[2], NormalCheck; 228 225 // get the two triangles 229 226 if (triangles.size() != 2) { 230 DoeLog(0) && (eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl);227 eLog() << Verbose(0) << "Baseline " << *this << " is connected to less than two triangles, Tesselation incomplete!" << endl; 231 228 return true; 232 229 } … … 236 233 BaseLineCenter.CopyVector(endpoints[0]->node->node); 237 234 BaseLineCenter.AddVector(endpoints[1]->node->node); 238 BaseLineCenter.Scale(1. /2.);235 BaseLineCenter.Scale(1./2.); 239 236 BaseLine.CopyVector(endpoints[0]->node->node); 240 237 BaseLine.SubtractVector(endpoints[1]->node->node); … … 244 241 NormalCheck.Zero(); 245 242 double sign = -1.; 246 int i =0;243 int i=0; 247 244 class BoundaryPointSet *node = NULL; 248 for (TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) {245 for(TriangleMap::const_iterator runner = triangles.begin(); runner != triangles.end(); runner++) { 249 246 //Log() << Verbose(0) << "INFO: NormalVector of " << *(runner->second) << " is " << runner->second->NormalVector << "." << endl; 250 247 NormalCheck.AddVector(&runner->second->NormalVector); … … 252 249 sign = -sign; 253 250 if (runner->second->NormalVector.NormSquared() > MYEPSILON) 254 BaseLineNormal.CopyVector(&runner->second->NormalVector); // yes, copy second on top of first251 BaseLineNormal.CopyVector(&runner->second->NormalVector); // yes, copy second on top of first 255 252 else { 256 DoeLog(0) && (eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl);253 eLog() << Verbose(0) << "Triangle " << *runner->second << " has zero normal vector!" << endl; 257 254 } 258 255 node = runner->second->GetThirdEndpoint(this); … … 261 258 helper[i].CopyVector(node->node->node); 262 259 helper[i].SubtractVector(&BaseLineCenter); 263 helper[i].MakeNormalVector(&BaseLine); // we want to compare the triangle's heights' angles!260 helper[i].MakeNormalVector(&BaseLine); // we want to compare the triangle's heights' angles! 264 261 //Log() << Verbose(0) << "INFO: Height vector with respect to baseline is " << helper[i] << "." << endl; 265 262 i++; 266 263 } else { 267 DoeLog(1) && (eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl);264 eLog() << Verbose(1) << "I cannot find third node in triangle, something's wrong." << endl; 268 265 return true; 269 266 } … … 271 268 //Log() << Verbose(0) << "INFO: BaselineNormal is " << BaseLineNormal << "." << endl; 272 269 if (NormalCheck.NormSquared() < MYEPSILON) { 273 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl);270 Log() << Verbose(0) << "ACCEPT: Normalvectors of both triangles are the same: convex." << endl; 274 271 return true; 275 272 } … … 277 274 double angle = GetAngle(helper[0], helper[1], BaseLineNormal); 278 275 if ((angle - M_PI) > -MYEPSILON) { 279 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl);276 Log() << Verbose(0) << "ACCEPT: Angle is greater than pi: convex." << endl; 280 277 return true; 281 278 } else { 282 DoLog(0) && (Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl);279 Log() << Verbose(0) << "REJECT: Angle is less than pi: concave." << endl; 283 280 return false; 284 281 } … … 291 288 bool BoundaryLineSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const 292 289 { 293 Info FunctionInfo(__func__);294 for (int i = 0; i < 2;i++)290 Info FunctionInfo(__func__); 291 for(int i=0;i<2;i++) 295 292 if (point == endpoints[i]) 296 293 return true; 297 294 return false; 298 } 299 ; 295 }; 300 296 301 297 /** Returns other endpoint of the line. … … 305 301 class BoundaryPointSet *BoundaryLineSet::GetOtherEndpoint(const BoundaryPointSet * const point) const 306 302 { 307 Info FunctionInfo(__func__);303 Info FunctionInfo(__func__); 308 304 if (endpoints[0] == point) 309 305 return endpoints[1]; … … 312 308 else 313 309 return NULL; 314 } 315 ; 310 }; 316 311 317 312 /** output operator for BoundaryLineSet. … … 319 314 * \param &a boundary line 320 315 */ 321 ostream & operator <<(ostream &ost, const BoundaryLineSet &a)316 ostream & operator <<(ostream &ost, const BoundaryLineSet &a) 322 317 { 323 318 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "]"; 324 319 return ost; 325 } 326 ; 320 }; 327 321 328 322 // ======================================== Triangles on Boundary ================================= … … 333 327 Nr(-1) 334 328 { 335 Info FunctionInfo(__func__);336 for (int i = 0; i < 3; i++) {337 endpoints[i] = NULL;338 lines[i] = NULL;339 }340 }341 ;329 Info FunctionInfo(__func__); 330 for (int i = 0; i < 3; i++) 331 { 332 endpoints[i] = NULL; 333 lines[i] = NULL; 334 } 335 }; 342 336 343 337 /** Constructor for BoundaryTriangleSet with three lines. … … 348 342 Nr(number) 349 343 { 350 Info FunctionInfo(__func__);344 Info FunctionInfo(__func__); 351 345 // set number 352 346 // set lines … … 360 354 // for all three lines 361 355 for (int j = 0; j < 2; j++) { // for both endpoints 362 OrderMap.insert(pair<int, class BoundaryPointSet *> (line[i]->endpoints[j]->Nr, line[i]->endpoints[j])); 356 OrderMap.insert(pair<int, class BoundaryPointSet *> ( 357 line[i]->endpoints[j]->Nr, line[i]->endpoints[j])); 363 358 // and we don't care whether insertion fails 364 359 } 365 360 // set endpoints 366 361 int Counter = 0; 367 DoLog(0) && (Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl);362 Log() << Verbose(0) << "New triangle " << Nr << " with end points: " << endl; 368 363 for (PointMap::iterator runner = OrderMap.begin(); runner != OrderMap.end(); runner++) { 369 364 endpoints[Counter] = runner->second; 370 DoLog(0) && (Log() << Verbose(0) << " " << *endpoints[Counter] << endl);365 Log() << Verbose(0) << " " << *endpoints[Counter] << endl; 371 366 Counter++; 372 367 } 373 368 if (Counter < 3) { 374 DoeLog(0) && (eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl);369 eLog() << Verbose(0) << "We have a triangle with only two distinct endpoints!" << endl; 375 370 performCriticalExit(); 376 371 } 377 } 378 ; 372 }; 379 373 380 374 /** Destructor of BoundaryTriangleSet. … … 384 378 BoundaryTriangleSet::~BoundaryTriangleSet() 385 379 { 386 Info FunctionInfo(__func__);380 Info FunctionInfo(__func__); 387 381 for (int i = 0; i < 3; i++) { 388 382 if (lines[i] != NULL) { … … 391 385 } 392 386 if (lines[i]->triangles.empty()) { 393 //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl;394 delete (lines[i]);395 lines[i] = NULL;387 //Log() << Verbose(0) << *lines[i] << " is no more attached to any triangle, erasing." << endl; 388 delete (lines[i]); 389 lines[i] = NULL; 396 390 } 397 391 } 398 392 } 399 393 //Log() << Verbose(0) << "Erasing triangle Nr." << Nr << " itself." << endl; 400 } 401 ; 394 }; 402 395 403 396 /** Calculates the normal vector for this triangle. … … 407 400 void BoundaryTriangleSet::GetNormalVector(const Vector &OtherVector) 408 401 { 409 Info FunctionInfo(__func__);402 Info FunctionInfo(__func__); 410 403 // get normal vector 411 404 NormalVector.MakeNormalVector(endpoints[0]->node->node, endpoints[1]->node->node, endpoints[2]->node->node); … … 414 407 if (NormalVector.ScalarProduct(&OtherVector) > 0.) 415 408 NormalVector.Scale(-1.); 416 DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl); 417 } 418 ; 409 Log() << Verbose(1) << "Normal Vector is " << NormalVector << "." << endl; 410 }; 419 411 420 412 /** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses. … … 437 429 438 430 if (!Intersection->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, MolCenter, x)) { 439 DoeLog(1) && (eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl);431 eLog() << Verbose(1) << "Alas! Intersection with plane failed - at least numerically - the intersection is not on the plane!" << endl; 440 432 return false; 441 433 } 442 434 443 DoLog(1) && (Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl);444 DoLog(1) && (Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl);445 DoLog(1) && (Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl);435 Log() << Verbose(1) << "INFO: Triangle is " << *this << "." << endl; 436 Log() << Verbose(1) << "INFO: Line is from " << *MolCenter << " to " << *x << "." << endl; 437 Log() << Verbose(1) << "INFO: Intersection is " << *Intersection << "." << endl; 446 438 447 439 if (Intersection->DistanceSquared(endpoints[0]->node->node) < MYEPSILON) { 448 DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl);440 Log() << Verbose(1) << "Intersection coindices with first endpoint." << endl; 449 441 return true; 450 } else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) {451 DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl);442 } else if (Intersection->DistanceSquared(endpoints[1]->node->node) < MYEPSILON) { 443 Log() << Verbose(1) << "Intersection coindices with second endpoint." << endl; 452 444 return true; 453 } else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) {454 DoLog(1) && (Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl);445 } else if (Intersection->DistanceSquared(endpoints[2]->node->node) < MYEPSILON) { 446 Log() << Verbose(1) << "Intersection coindices with third endpoint." << endl; 455 447 return true; 456 448 } 457 449 // Calculate cross point between one baseline and the line from the third endpoint to intersection 458 int i =0;450 int i=0; 459 451 do { 460 if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i % 3]->node->node, endpoints[(i + 1) % 3]->node->node, endpoints[(i + 2) %3]->node->node, Intersection, &NormalVector)) {461 helper.CopyVector(endpoints[(i + 1) %3]->node->node);462 helper.SubtractVector(endpoints[i %3]->node->node);463 CrossPoint.SubtractVector(endpoints[i % 3]->node->node);// cross point was returned as absolute vector464 const double s = CrossPoint.ScalarProduct(&helper) /helper.NormSquared();465 DoLog(1) && (Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl);466 if ((s < -MYEPSILON) || ((s -1.) > MYEPSILON)) {467 DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl);468 i =4;452 if (CrossPoint.GetIntersectionOfTwoLinesOnPlane(endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node, endpoints[(i+2)%3]->node->node, Intersection, &NormalVector)) { 453 helper.CopyVector(endpoints[(i+1)%3]->node->node); 454 helper.SubtractVector(endpoints[i%3]->node->node); 455 CrossPoint.SubtractVector(endpoints[i%3]->node->node); // cross point was returned as absolute vector 456 const double s = CrossPoint.ScalarProduct(&helper)/helper.NormSquared(); 457 Log() << Verbose(1) << "INFO: Factor s is " << s << "." << endl; 458 if ((s < -MYEPSILON) || ((s-1.) > MYEPSILON)) { 459 Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << "outside of triangle." << endl; 460 i=4; 469 461 break; 470 462 } 471 463 i++; 472 } else 464 } else 473 465 break; 474 } while (i <3);475 if (i ==3) {476 DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl);466 } while (i<3); 467 if (i==3) { 468 Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " inside of triangle." << endl; 477 469 return true; 478 470 } else { 479 DoLog(1) && (Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl);471 Log() << Verbose(1) << "INFO: Crosspoint " << CrossPoint << " outside of triangle." << endl; 480 472 return false; 481 473 } 482 } 483 ; 484 485 /** Finds the point on the triangle to the point \a *x. 486 * We call Vector::GetIntersectionWithPlane() with \a * and the center of the triangle to receive an intersection point. 487 * Then we check the in-plane part (the part projected down onto plane). We check whether it crosses one of the 488 * boundary lines. If it does, we return this intersection as closest point, otherwise the projected point down. 474 }; 475 476 /** Finds the point on the triangle \a *BTS through which the line defined by \a *MolCenter and \a *x crosses. 477 * We call Vector::GetIntersectionWithPlane() to receive the intersection point with the plane 489 478 * Thus we test if it's really on the plane and whether it's inside the triangle on the plane or not. 490 479 * The latter is done as follows: We calculate the cross point of one of the triangle's baseline with the line … … 501 490 502 491 // 1. get intersection with plane 503 DoLog(1) && (Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl);492 Log() << Verbose(1) << "INFO: Looking for closest point of triangle " << *this << " to " << *x << "." << endl; 504 493 GetCenter(&Direction); 505 494 if (!ClosestPoint->GetIntersectionWithPlane(&NormalVector, endpoints[0]->node->node, x, &Direction)) { … … 510 499 Vector InPlane; 511 500 InPlane.CopyVector(x); 512 InPlane.SubtractVector(ClosestPoint); // points from plane intersection to straight-down point501 InPlane.SubtractVector(ClosestPoint); // points from plane intersection to straight-down point 513 502 InPlane.ProjectOntoPlane(&NormalVector); 514 503 InPlane.AddVector(ClosestPoint); 515 504 516 DoLog(2) && (Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl);517 DoLog(2) && (Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl);518 DoLog(2) && (Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl);505 Log() << Verbose(2) << "INFO: Triangle is " << *this << "." << endl; 506 Log() << Verbose(2) << "INFO: Line is from " << Direction << " to " << *x << "." << endl; 507 Log() << Verbose(2) << "INFO: In-plane part is " << InPlane << "." << endl; 519 508 520 509 // Calculate cross point between one baseline and the desired point such that distance is shortest … … 524 513 Vector CrossPoint[3]; 525 514 Vector helper; 526 for (int i = 0; i < 3;i++) {515 for (int i=0;i<3;i++) { 527 516 // treat direction of line as normal of a (cut)plane and the desired point x as the plane offset, the intersect line with point 528 Direction.CopyVector(endpoints[(i + 1) %3]->node->node);529 Direction.SubtractVector(endpoints[i %3]->node->node);517 Direction.CopyVector(endpoints[(i+1)%3]->node->node); 518 Direction.SubtractVector(endpoints[i%3]->node->node); 530 519 // calculate intersection, line can never be parallel to Direction (is the same vector as PlaneNormal); 531 CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i % 3]->node->node, endpoints[(i + 1) %3]->node->node);520 CrossPoint[i].GetIntersectionWithPlane(&Direction, &InPlane, endpoints[i%3]->node->node, endpoints[(i+1)%3]->node->node); 532 521 CrossDirection[i].CopyVector(&CrossPoint[i]); 533 522 CrossDirection[i].SubtractVector(&InPlane); 534 CrossPoint[i].SubtractVector(endpoints[i % 3]->node->node);// cross point was returned as absolute vector535 const double s = CrossPoint[i].ScalarProduct(&Direction) /Direction.NormSquared();536 DoLog(2) && (Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl);537 if ((s >= -MYEPSILON) && ((s -1.) <= MYEPSILON)) {538 CrossPoint[i].AddVector(endpoints[i % 3]->node->node);// make cross point absolute again539 DoLog(2) && (Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i % 3]->node->node << " and " << *endpoints[(i + 1) % 3]->node->node << "." << endl);523 CrossPoint[i].SubtractVector(endpoints[i%3]->node->node); // cross point was returned as absolute vector 524 const double s = CrossPoint[i].ScalarProduct(&Direction)/Direction.NormSquared(); 525 Log() << Verbose(2) << "INFO: Factor s is " << s << "." << endl; 526 if ((s >= -MYEPSILON) && ((s-1.) <= MYEPSILON)) { 527 CrossPoint[i].AddVector(endpoints[i%3]->node->node); // make cross point absolute again 528 Log() << Verbose(2) << "INFO: Crosspoint is " << CrossPoint[i] << ", intersecting BoundaryLine between " << *endpoints[i%3]->node->node << " and " << *endpoints[(i+1)%3]->node->node << "." << endl; 540 529 const double distance = CrossPoint[i].DistanceSquared(x); 541 530 if ((ShortestDistance < 0.) || (ShortestDistance > distance)) { … … 547 536 } 548 537 InsideFlag = true; 549 for (int i = 0; i < 3; i++) { 550 const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i + 1) % 3]); 551 const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i + 2) % 3]); 552 ; 553 if ((sign > -MYEPSILON) && (othersign > -MYEPSILON)) // have different sign 538 for (int i=0;i<3;i++) { 539 const double sign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+1)%3]); 540 const double othersign = CrossDirection[i].ScalarProduct(&CrossDirection[(i+2)%3]);; 541 if ((sign > -MYEPSILON) && (othersign > -MYEPSILON)) // have different sign 554 542 InsideFlag = false; 555 543 } … … 557 545 ClosestPoint->CopyVector(&InPlane); 558 546 ShortestDistance = InPlane.DistanceSquared(x); 559 } else { // also check endnodes560 for (int i = 0; i < 3;i++) {547 } else { // also check endnodes 548 for (int i=0;i<3;i++) { 561 549 const double distance = x->DistanceSquared(endpoints[i]->node->node); 562 550 if ((ShortestDistance < 0.) || (ShortestDistance > distance)) { … … 566 554 } 567 555 } 568 DoLog(1) && (Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl);556 Log() << Verbose(1) << "INFO: Closest Point is " << *ClosestPoint << " with shortest squared distance is " << ShortestDistance << "." << endl; 569 557 return ShortestDistance; 570 } 571 ; 558 }; 572 559 573 560 /** Checks whether lines is any of the three boundary lines this triangle contains. … … 577 564 bool BoundaryTriangleSet::ContainsBoundaryLine(const BoundaryLineSet * const line) const 578 565 { 579 Info FunctionInfo(__func__);580 for (int i = 0; i < 3;i++)566 Info FunctionInfo(__func__); 567 for(int i=0;i<3;i++) 581 568 if (line == lines[i]) 582 569 return true; 583 570 return false; 584 } 585 ; 571 }; 586 572 587 573 /** Checks whether point is any of the three endpoints this triangle contains. … … 591 577 bool BoundaryTriangleSet::ContainsBoundaryPoint(const BoundaryPointSet * const point) const 592 578 { 593 Info FunctionInfo(__func__);594 for (int i = 0; i < 3;i++)579 Info FunctionInfo(__func__); 580 for(int i=0;i<3;i++) 595 581 if (point == endpoints[i]) 596 582 return true; 597 583 return false; 598 } 599 ; 584 }; 600 585 601 586 /** Checks whether point is any of the three endpoints this triangle contains. … … 605 590 bool BoundaryTriangleSet::ContainsBoundaryPoint(const TesselPoint * const point) const 606 591 { 607 Info FunctionInfo(__func__);608 for (int i = 0; i < 3;i++)592 Info FunctionInfo(__func__); 593 for(int i=0;i<3;i++) 609 594 if (point == endpoints[i]->node) 610 595 return true; 611 596 return false; 612 } 613 ; 597 }; 614 598 615 599 /** Checks whether three given \a *Points coincide with triangle's endpoints. … … 619 603 bool BoundaryTriangleSet::IsPresentTupel(const BoundaryPointSet * const Points[3]) const 620 604 { 621 Info FunctionInfo(__func__); 622 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking " << Points[0] << "," << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl); 623 return (((endpoints[0] == Points[0]) || (endpoints[0] == Points[1]) || (endpoints[0] == Points[2])) && ((endpoints[1] == Points[0]) || (endpoints[1] == Points[1]) || (endpoints[1] == Points[2])) && ((endpoints[2] == Points[0]) || (endpoints[2] == Points[1]) || (endpoints[2] == Points[2]) 624 625 )); 626 } 627 ; 605 Info FunctionInfo(__func__); 606 Log() << Verbose(1) << "INFO: Checking " << Points[0] << "," << Points[1] << "," << Points[2] << " against " << endpoints[0] << "," << endpoints[1] << "," << endpoints[2] << "." << endl; 607 return (((endpoints[0] == Points[0]) 608 || (endpoints[0] == Points[1]) 609 || (endpoints[0] == Points[2]) 610 ) && ( 611 (endpoints[1] == Points[0]) 612 || (endpoints[1] == Points[1]) 613 || (endpoints[1] == Points[2]) 614 ) && ( 615 (endpoints[2] == Points[0]) 616 || (endpoints[2] == Points[1]) 617 || (endpoints[2] == Points[2]) 618 619 )); 620 }; 628 621 629 622 /** Checks whether three given \a *Points coincide with triangle's endpoints. … … 633 626 bool BoundaryTriangleSet::IsPresentTupel(const BoundaryTriangleSet * const T) const 634 627 { 635 Info FunctionInfo(__func__); 636 return (((endpoints[0] == T->endpoints[0]) || (endpoints[0] == T->endpoints[1]) || (endpoints[0] == T->endpoints[2])) && ((endpoints[1] == T->endpoints[0]) || (endpoints[1] == T->endpoints[1]) || (endpoints[1] == T->endpoints[2])) && ((endpoints[2] == T->endpoints[0]) || (endpoints[2] == T->endpoints[1]) || (endpoints[2] == T->endpoints[2]) 637 638 )); 639 } 640 ; 628 Info FunctionInfo(__func__); 629 return (((endpoints[0] == T->endpoints[0]) 630 || (endpoints[0] == T->endpoints[1]) 631 || (endpoints[0] == T->endpoints[2]) 632 ) && ( 633 (endpoints[1] == T->endpoints[0]) 634 || (endpoints[1] == T->endpoints[1]) 635 || (endpoints[1] == T->endpoints[2]) 636 ) && ( 637 (endpoints[2] == T->endpoints[0]) 638 || (endpoints[2] == T->endpoints[1]) 639 || (endpoints[2] == T->endpoints[2]) 640 641 )); 642 }; 641 643 642 644 /** Returns the endpoint which is not contained in the given \a *line. … … 646 648 class BoundaryPointSet *BoundaryTriangleSet::GetThirdEndpoint(const BoundaryLineSet * const line) const 647 649 { 648 Info FunctionInfo(__func__);650 Info FunctionInfo(__func__); 649 651 // sanity check 650 652 if (!ContainsBoundaryLine(line)) 651 653 return NULL; 652 for (int i = 0; i < 3;i++)654 for(int i=0;i<3;i++) 653 655 if (!line->ContainsBoundaryPoint(endpoints[i])) 654 656 return endpoints[i]; 655 657 // actually, that' impossible :) 656 658 return NULL; 657 } 658 ; 659 }; 659 660 660 661 /** Calculates the center point of the triangle. … … 664 665 void BoundaryTriangleSet::GetCenter(Vector * const center) const 665 666 { 666 Info FunctionInfo(__func__);667 Info FunctionInfo(__func__); 667 668 center->Zero(); 668 for (int i = 0; i < 3;i++)669 for(int i=0;i<3;i++) 669 670 center->AddVector(endpoints[i]->node->node); 670 center->Scale(1. /3.);671 DoLog(1) && (Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl);671 center->Scale(1./3.); 672 Log() << Verbose(1) << "INFO: Center is at " << *center << "." << endl; 672 673 } 673 674 … … 679 680 { 680 681 ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << "," << a.endpoints[1]->node->Name << "," << a.endpoints[2]->node->Name << "]"; 681 // ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << ","682 // << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]";682 // ost << "[" << a.Nr << "|" << a.endpoints[0]->node->Name << " at " << *a.endpoints[0]->node->node << "," 683 // << a.endpoints[1]->node->Name << " at " << *a.endpoints[1]->node->node << "," << a.endpoints[2]->node->Name << " at " << *a.endpoints[2]->node->node << "]"; 683 684 return ost; 684 } 685 ; 685 }; 686 686 687 687 // ======================================== Polygons on Boundary ================================= … … 693 693 { 694 694 Info FunctionInfo(__func__); 695 } 696 ; 695 }; 697 696 698 697 /** Destructor of BoundaryPolygonSet. … … 704 703 Info FunctionInfo(__func__); 705 704 endpoints.clear(); 706 DoLog(1) && (Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl); 707 } 708 ; 705 Log() << Verbose(1) << "Erasing polygon Nr." << Nr << " itself." << endl; 706 }; 709 707 710 708 /** Calculates the normal vector for this triangle. … … 720 718 Vector *TotalNormal = new Vector; 721 719 PointSet::const_iterator Runner[3]; 722 for (int i = 0; i <3; i++) {720 for (int i=0;i<3; i++) { 723 721 Runner[i] = endpoints.begin(); 724 for (int j = 0; j <i; j++) { // go as much further722 for (int j = 0; j<i; j++) { // go as much further 725 723 Runner[i]++; 726 724 if (Runner[i] == endpoints.end()) { 727 DoeLog(0) && (eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl);725 eLog() << Verbose(0) << "There are less than three endpoints in the polygon!" << endl; 728 726 performCriticalExit(); 729 727 } … … 731 729 } 732 730 TotalNormal->Zero(); 733 int counter =0;734 for (; Runner[2] != endpoints.end(); ) {731 int counter=0; 732 for (; Runner[2] != endpoints.end(); ) { 735 733 TemporaryNormal.MakeNormalVector((*Runner[0])->node->node, (*Runner[1])->node->node, (*Runner[2])->node->node); 736 for (int i = 0; i < 3;i++) // increase each of them734 for (int i=0;i<3;i++) // increase each of them 737 735 Runner[i]++; 738 736 TotalNormal->AddVector(&TemporaryNormal); 739 737 } 740 TotalNormal->Scale(1. / (double)counter);738 TotalNormal->Scale(1./(double)counter); 741 739 742 740 // make it always point inward (any offset vector onto plane projected onto normal vector suffices) 743 741 if (TotalNormal->ScalarProduct(&OtherVector) > 0.) 744 742 TotalNormal->Scale(-1.); 745 DoLog(1) && (Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl);743 Log() << Verbose(1) << "Normal Vector is " << *TotalNormal << "." << endl; 746 744 747 745 return TotalNormal; 748 } 749 ; 746 }; 750 747 751 748 /** Calculates the center point of the triangle. … … 758 755 center->Zero(); 759 756 int counter = 0; 760 for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {757 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) { 761 758 center->AddVector((*Runner)->node->node); 762 759 counter++; 763 760 } 764 center->Scale(1. / (double)counter);765 DoLog(1) && (Log() << Verbose(1) << "Center is at " << *center << "." << endl);761 center->Scale(1./(double)counter); 762 Log() << Verbose(1) << "Center is at " << *center << "." << endl; 766 763 } 767 764 … … 774 771 Info FunctionInfo(__func__); 775 772 return ContainsPresentTupel(triangle->endpoints, 3); 776 } 777 ; 773 }; 778 774 779 775 /** Checks whether the polygons contains both endpoints of the line. … … 785 781 Info FunctionInfo(__func__); 786 782 return ContainsPresentTupel(line->endpoints, 2); 787 } 788 ; 783 }; 789 784 790 785 /** Checks whether point is any of the three endpoints this triangle contains. … … 795 790 { 796 791 Info FunctionInfo(__func__); 797 for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {798 DoLog(0) && (Log() << Verbose(0) << "Checking against " << **Runner << endl);792 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) { 793 Log() << Verbose(0) << "Checking against " << **Runner << endl; 799 794 if (point == (*Runner)) { 800 DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);795 Log() << Verbose(0) << " Contained." << endl; 801 796 return true; 802 797 } 803 798 } 804 DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);799 Log() << Verbose(0) << " Not contained." << endl; 805 800 return false; 806 } 807 ; 801 }; 808 802 809 803 /** Checks whether point is any of the three endpoints this triangle contains. … … 814 808 { 815 809 Info FunctionInfo(__func__); 816 for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)810 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) 817 811 if (point == (*Runner)->node) { 818 DoLog(0) && (Log() << Verbose(0) << " Contained." << endl);812 Log() << Verbose(0) << " Contained." << endl; 819 813 return true; 820 814 } 821 DoLog(0) && (Log() << Verbose(0) << " Not contained." << endl);815 Log() << Verbose(0) << " Not contained." << endl; 822 816 return false; 823 } 824 ; 817 }; 825 818 826 819 /** Checks whether given array of \a *Points coincide with polygons's endpoints. … … 833 826 Info FunctionInfo(__func__); 834 827 int counter = 0; 835 DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);836 for (int i = 0; i < dim;i++) {837 DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl);828 Log() << Verbose(1) << "Polygon is " << *this << endl; 829 for(int i=0;i<dim;i++) { 830 Log() << Verbose(1) << " Testing endpoint " << *Points[i] << endl; 838 831 if (ContainsBoundaryPoint(Points[i])) { 839 832 counter++; … … 845 838 else 846 839 return false; 847 } 848 ; 840 }; 849 841 850 842 /** Checks whether given PointList coincide with polygons's endpoints. … … 856 848 Info FunctionInfo(__func__); 857 849 size_t counter = 0; 858 DoLog(1) && (Log() << Verbose(1) << "Polygon is " << *this << endl);859 for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) {860 DoLog(1) && (Log() << Verbose(1) << " Testing endpoint " << **Runner << endl);850 Log() << Verbose(1) << "Polygon is " << *this << endl; 851 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) { 852 Log() << Verbose(1) << " Testing endpoint " << **Runner << endl; 861 853 if (ContainsBoundaryPoint(*Runner)) 862 854 counter++; … … 867 859 else 868 860 return false; 869 } 870 ; 861 }; 871 862 872 863 /** Checks whether given set of \a *Points coincide with polygons's endpoints. … … 876 867 bool BoundaryPolygonSet::ContainsPresentTupel(const BoundaryPolygonSet * const P) const 877 868 { 878 return ContainsPresentTupel((const PointSet) P->endpoints); 879 } 880 ; 869 return ContainsPresentTupel((const PointSet)P->endpoints); 870 }; 881 871 882 872 /** Gathers all the endpoints' triangles in a unique set. … … 886 876 { 887 877 Info FunctionInfo(__func__); 888 pair <TriangleSet::iterator, bool> Tester;878 pair <TriangleSet::iterator, bool> Tester; 889 879 TriangleSet *triangles = new TriangleSet; 890 880 891 for (PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++)892 for (LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++)893 for (TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) {881 for(PointSet::const_iterator Runner = endpoints.begin(); Runner != endpoints.end(); Runner++) 882 for(LineMap::const_iterator Walker = (*Runner)->lines.begin(); Walker != (*Runner)->lines.end(); Walker++) 883 for(TriangleMap::const_iterator Sprinter = (Walker->second)->triangles.begin(); Sprinter != (Walker->second)->triangles.end(); Sprinter++) { 894 884 //Log() << Verbose(0) << " Testing triangle " << *(Sprinter->second) << endl; 895 885 if (ContainsBoundaryTriangle(Sprinter->second)) { 896 886 Tester = triangles->insert(Sprinter->second); 897 887 if (Tester.second) 898 DoLog(0) && (Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl);888 Log() << Verbose(0) << "Adding triangle " << *(Sprinter->second) << endl; 899 889 } 900 890 } 901 891 902 DoLog(1) && (Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl);892 Log() << Verbose(1) << "The Polygon of " << endpoints.size() << " endpoints has " << triangles->size() << " unique triangles in total." << endl; 903 893 return triangles; 904 } 905 ; 894 }; 906 895 907 896 /** Fills the endpoints of this polygon from the triangles attached to \a *line. … … 912 901 { 913 902 Info FunctionInfo(__func__); 914 pair <PointSet::iterator, bool> Tester;903 pair <PointSet::iterator, bool> Tester; 915 904 if (line == NULL) 916 905 return false; 917 DoLog(1) && (Log() << Verbose(1) << "Filling polygon from line " << *line << endl);918 for (TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) {919 for (int i = 0; i < 3;i++) {906 Log() << Verbose(1) << "Filling polygon from line " << *line << endl; 907 for(TriangleMap::const_iterator Runner = line->triangles.begin(); Runner != line->triangles.end(); Runner++) { 908 for (int i=0;i<3;i++) { 920 909 Tester = endpoints.insert((Runner->second)->endpoints[i]); 921 910 if (Tester.second) 922 DoLog(1) && (Log() << Verbose(1) << " Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl);911 Log() << Verbose(1) << " Inserting endpoint " << *((Runner->second)->endpoints[i]) << endl; 923 912 } 924 913 } 925 914 926 915 return true; 927 } 928 ; 916 }; 929 917 930 918 /** output operator for BoundaryPolygonSet. … … 935 923 { 936 924 ost << "[" << a.Nr << "|"; 937 for (PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) {938 ost << (*Runner)->node->Name;939 Runner++;940 if (Runner != a.endpoints.end())941 ost << ",";942 } 943 ost << "]";925 for(PointSet::const_iterator Runner = a.endpoints.begin(); Runner != a.endpoints.end();) { 926 ost << (*Runner)->node->Name; 927 Runner++; 928 if (Runner != a.endpoints.end()) 929 ost << ","; 930 } 931 ost<< "]"; 944 932 return ost; 945 } 946 ; 933 }; 947 934 948 935 // =========================================================== class TESSELPOINT =========================================== … … 955 942 node = NULL; 956 943 nr = -1; 957 Name = NULL; 958 } 959 ; 944 Name = NULL; 945 }; 960 946 961 947 /** Destructor for class TesselPoint. … … 964 950 { 965 951 //Info FunctionInfo(__func__); 966 } 967 ; 952 }; 968 953 969 954 /** Prints LCNode to screen. 970 955 */ 971 ostream & operator << (ostream &ost, const TesselPoint &a)956 ostream & operator << (ostream &ost, const TesselPoint &a) 972 957 { 973 958 ost << "[" << (a.Name) << "|" << a.Name << " at " << *a.node << "]"; 974 959 return ost; 975 } 976 ; 960 }; 977 961 978 962 /** Prints LCNode to screen. 979 963 */ 980 ostream & TesselPoint::operator << (ostream &ost)981 { 982 Info FunctionInfo(__func__);964 ostream & TesselPoint::operator << (ostream &ost) 965 { 966 Info FunctionInfo(__func__); 983 967 ost << "[" << (nr) << "|" << this << "]"; 984 968 return ost; 985 } 986 ; 969 }; 970 987 971 988 972 // =========================================================== class POINTCLOUD ============================================ … … 992 976 PointCloud::PointCloud() 993 977 { 994 //Info FunctionInfo(__func__); 995 } 996 ; 978 //Info FunctionInfo(__func__); 979 }; 997 980 998 981 /** Destructor for class PointCloud. … … 1000 983 PointCloud::~PointCloud() 1001 984 { 1002 //Info FunctionInfo(__func__); 1003 } 1004 ; 985 //Info FunctionInfo(__func__); 986 }; 1005 987 1006 988 // ============================ CandidateForTesselation ============================= … … 1008 990 /** Constructor of class CandidateForTesselation. 1009 991 */ 1010 CandidateForTesselation::CandidateForTesselation(BoundaryLineSet* line) : 1011 BaseLine(line), ThirdPoint(NULL), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI) 1012 { 1013 Info FunctionInfo(__func__); 1014 } 1015 ; 992 CandidateForTesselation::CandidateForTesselation (BoundaryLineSet* line) : 993 BaseLine(line), 994 ShortestAngle(2.*M_PI), 995 OtherShortestAngle(2.*M_PI) 996 { 997 Info FunctionInfo(__func__); 998 }; 999 1016 1000 1017 1001 /** Constructor of class CandidateForTesselation. 1018 1002 */ 1019 CandidateForTesselation::CandidateForTesselation(TesselPoint *candidate, BoundaryLineSet* line, BoundaryPointSet* point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) : 1020 BaseLine(line), ThirdPoint(point), T(NULL), ShortestAngle(2. * M_PI), OtherShortestAngle(2. * M_PI) 1021 { 1022 Info FunctionInfo(__func__); 1003 CandidateForTesselation::CandidateForTesselation (TesselPoint *candidate, BoundaryLineSet* line, Vector OptCandidateCenter, Vector OtherOptCandidateCenter) : 1004 BaseLine(line), 1005 ShortestAngle(2.*M_PI), 1006 OtherShortestAngle(2.*M_PI) 1007 { 1008 Info FunctionInfo(__func__); 1023 1009 OptCenter.CopyVector(&OptCandidateCenter); 1024 1010 OtherOptCenter.CopyVector(&OtherOptCandidateCenter); 1025 } 1026 ; 1011 }; 1027 1012 1028 1013 /** Destructor for class CandidateForTesselation. 1029 1014 */ 1030 CandidateForTesselation::~CandidateForTesselation() 1031 { 1032 } 1033 ; 1034 1035 /** Checks validity of a given sphere of a candidate line. 1036 * Sphere must touch all candidates and the baseline endpoints and there must be no other atoms inside. 1037 * \param RADIUS radius of sphere 1038 * \param *LC LinkedCell structure with other atoms 1039 * \return true - sphere is valid, false - sphere contains other points 1040 */ 1041 bool CandidateForTesselation::CheckValidity(const double RADIUS, const LinkedCell *LC) const 1042 { 1043 Info FunctionInfo(__func__); 1044 1045 const double radiusSquared = RADIUS * RADIUS; 1046 list<const Vector *> VectorList; 1047 VectorList.push_back(&OptCenter); 1048 //VectorList.push_back(&OtherOptCenter); // don't check the other (wrong) center 1049 1050 if (!pointlist.empty()) 1051 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains candidate list and baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl); 1052 else 1053 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere with no candidates contains baseline " << *BaseLine->endpoints[0] << "<->" << *BaseLine->endpoints[1] << " only ..." << endl); 1054 // check baseline for OptCenter and OtherOptCenter being on sphere's surface 1055 for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { 1056 for (int i = 0; i < 2; i++) { 1057 const double distance = fabs((*VRunner)->DistanceSquared(BaseLine->endpoints[i]->node->node) - radiusSquared); 1058 if (distance > HULLEPSILON) { 1059 DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << *BaseLine->endpoints[i] << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl); 1060 return false; 1061 } 1062 } 1063 } 1064 1065 // check Candidates for OptCenter and OtherOptCenter being on sphere's surface 1066 for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) { 1067 const TesselPoint *Walker = *Runner; 1068 for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { 1069 const double distance = fabs((*VRunner)->DistanceSquared(Walker->node) - radiusSquared); 1070 if (distance > HULLEPSILON) { 1071 DoeLog(1) && (eLog() << Verbose(1) << "Candidate " << *Walker << " is out of sphere at " << *(*VRunner) << " by " << distance << "." << endl); 1072 return false; 1073 } else { 1074 DoLog(1) && (Log() << Verbose(1) << "Candidate " << *Walker << " is inside by " << distance << "." << endl); 1075 } 1076 } 1077 } 1078 1079 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl); 1080 bool flag = true; 1081 for (list<const Vector *>::const_iterator VRunner = VectorList.begin(); VRunner != VectorList.end(); ++VRunner) { 1082 // get all points inside the sphere 1083 TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, (*VRunner)); 1084 1085 DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << OtherOptCenter << ":" << endl); 1086 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) 1087 DoLog(1) && (Log() << Verbose(1) << " " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&OtherOptCenter) << "." << endl); 1088 1089 // remove baseline's endpoints and candidates 1090 for (int i = 0; i < 2; i++) { 1091 DoLog(1) && (Log() << Verbose(1) << "INFO: removing baseline tesselpoint " << *BaseLine->endpoints[i]->node << "." << endl); 1092 ListofPoints->remove(BaseLine->endpoints[i]->node); 1093 } 1094 for (TesselPointList::const_iterator Runner = pointlist.begin(); Runner != pointlist.end(); ++Runner) { 1095 DoLog(1) && (Log() << Verbose(1) << "INFO: removing candidate tesselpoint " << *(*Runner) << "." << endl); 1096 ListofPoints->remove(*Runner); 1097 } 1098 if (!ListofPoints->empty()) { 1099 DoeLog(1) && (eLog() << Verbose(1) << "CheckValidity: There are still " << ListofPoints->size() << " points inside the sphere." << endl); 1100 flag = false; 1101 DoeLog(1) && (eLog() << Verbose(1) << "External atoms inside of sphere at " << *(*VRunner) << ":" << endl); 1102 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) 1103 DoeLog(1) && (eLog() << Verbose(1) << " " << *(*Runner) << endl); 1104 } 1105 delete (ListofPoints); 1106 1107 // check with animate_sphere.tcl VMD script 1108 if (ThirdPoint != NULL) { 1109 DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " " << ThirdPoint->Nr + 1 << " " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl); 1110 } else { 1111 DoLog(1) && (Log() << Verbose(1) << "Check by: ... missing third point ..." << endl); 1112 DoLog(1) && (Log() << Verbose(1) << "Check by: animate_sphere 0 " << BaseLine->endpoints[0]->Nr + 1 << " " << BaseLine->endpoints[1]->Nr + 1 << " ??? " << RADIUS << " " << OldCenter.x[0] << " " << OldCenter.x[1] << " " << OldCenter.x[2] << " " << (*VRunner)->x[0] << " " << (*VRunner)->x[1] << " " << (*VRunner)->x[2] << endl); 1113 } 1114 } 1115 return flag; 1116 } 1117 ; 1015 CandidateForTesselation::~CandidateForTesselation() { 1016 BaseLine = NULL; 1017 }; 1118 1018 1119 1019 /** output operator for CandidateForTesselation. … … 1121 1021 * \param &a boundary line 1122 1022 */ 1123 ostream & operator <<(ostream &ost, const CandidateForTesselation &a)1023 ostream & operator <<(ostream &ost, const CandidateForTesselation &a) 1124 1024 { 1125 1025 ost << "[" << a.BaseLine->Nr << "|" << a.BaseLine->endpoints[0]->node->Name << "," << a.BaseLine->endpoints[1]->node->Name << "] with "; … … 1134 1034 for (TesselPointList::const_iterator Runner = a.pointlist.begin(); Runner != a.pointlist.end(); Runner++) 1135 1035 ost << *(*Runner) << " "; 1136 ost << " at angle " << (a.ShortestAngle) << ".";1036 ost << " at angle " << (a.ShortestAngle)<< "."; 1137 1037 } 1138 1038 1139 1039 return ost; 1040 }; 1041 1042 1043 // =========================================================== class TESSELATION =========================================== 1044 1045 /** Constructor of class Tesselation. 1046 */ 1047 Tesselation::Tesselation() : 1048 PointsOnBoundaryCount(0), 1049 LinesOnBoundaryCount(0), 1050 TrianglesOnBoundaryCount(0), 1051 LastTriangle(NULL), 1052 TriangleFilesWritten(0), 1053 InternalPointer(PointsOnBoundary.begin()) 1054 { 1055 Info FunctionInfo(__func__); 1140 1056 } 1141 1057 ; 1142 1058 1143 // =========================================================== class TESSELATION ===========================================1144 1145 /** Constructor of class Tesselation.1146 */1147 Tesselation::Tesselation() :1148 PointsOnBoundaryCount(0), LinesOnBoundaryCount(0), TrianglesOnBoundaryCount(0), LastTriangle(NULL), TriangleFilesWritten(0), InternalPointer(PointsOnBoundary.begin())1149 {1150 Info FunctionInfo(__func__);1151 }1152 ;1153 1154 1059 /** Destructor of class Tesselation. 1155 1060 * We have to free all points, lines and triangles. … … 1157 1062 Tesselation::~Tesselation() 1158 1063 { 1159 Info FunctionInfo(__func__);1160 DoLog(0) && (Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl);1064 Info FunctionInfo(__func__); 1065 Log() << Verbose(0) << "Free'ing TesselStruct ... " << endl; 1161 1066 for (TriangleMap::iterator runner = TrianglesOnBoundary.begin(); runner != TrianglesOnBoundary.end(); runner++) { 1162 1067 if (runner->second != NULL) { … … 1164 1069 runner->second = NULL; 1165 1070 } else 1166 DoeLog(1) && (eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl);1167 } 1168 DoLog(0) && (Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl);1071 eLog() << Verbose(1) << "The triangle " << runner->first << " has already been free'd." << endl; 1072 } 1073 Log() << Verbose(0) << "This envelope was written to file " << TriangleFilesWritten << " times(s)." << endl; 1169 1074 } 1170 1075 ; … … 1172 1077 /** PointCloud implementation of GetCenter 1173 1078 * Uses PointsOnBoundary and STL stuff. 1174 */ 1079 */ 1175 1080 Vector * Tesselation::GetCenter(ofstream *out) const 1176 1081 { 1177 Info FunctionInfo(__func__);1178 Vector *Center = new Vector(0., 0.,0.);1179 int num =0;1082 Info FunctionInfo(__func__); 1083 Vector *Center = new Vector(0.,0.,0.); 1084 int num=0; 1180 1085 for (GoToFirst(); (!IsEnd()); GoToNext()) { 1181 1086 Center->AddVector(GetPoint()->node); 1182 1087 num++; 1183 1088 } 1184 Center->Scale(1. /num);1089 Center->Scale(1./num); 1185 1090 return Center; 1186 } 1187 ; 1091 }; 1188 1092 1189 1093 /** PointCloud implementation of GoPoint 1190 1094 * Uses PointsOnBoundary and STL stuff. 1191 */ 1095 */ 1192 1096 TesselPoint * Tesselation::GetPoint() const 1193 1097 { 1194 Info FunctionInfo(__func__);1098 Info FunctionInfo(__func__); 1195 1099 return (InternalPointer->second->node); 1196 } 1197 ; 1100 }; 1198 1101 1199 1102 /** PointCloud implementation of GetTerminalPoint. 1200 1103 * Uses PointsOnBoundary and STL stuff. 1201 */ 1104 */ 1202 1105 TesselPoint * Tesselation::GetTerminalPoint() const 1203 1106 { 1204 Info FunctionInfo(__func__);1107 Info FunctionInfo(__func__); 1205 1108 PointMap::const_iterator Runner = PointsOnBoundary.end(); 1206 1109 Runner--; 1207 1110 return (Runner->second->node); 1208 } 1209 ; 1111 }; 1210 1112 1211 1113 /** PointCloud implementation of GoToNext. 1212 1114 * Uses PointsOnBoundary and STL stuff. 1213 */ 1115 */ 1214 1116 void Tesselation::GoToNext() const 1215 1117 { 1216 Info FunctionInfo(__func__);1118 Info FunctionInfo(__func__); 1217 1119 if (InternalPointer != PointsOnBoundary.end()) 1218 1120 InternalPointer++; 1219 } 1220 ; 1121 }; 1221 1122 1222 1123 /** PointCloud implementation of GoToPrevious. 1223 1124 * Uses PointsOnBoundary and STL stuff. 1224 */ 1125 */ 1225 1126 void Tesselation::GoToPrevious() const 1226 1127 { 1227 Info FunctionInfo(__func__);1128 Info FunctionInfo(__func__); 1228 1129 if (InternalPointer != PointsOnBoundary.begin()) 1229 1130 InternalPointer--; 1230 } 1231 ; 1131 }; 1232 1132 1233 1133 /** PointCloud implementation of GoToFirst. 1234 1134 * Uses PointsOnBoundary and STL stuff. 1235 */ 1135 */ 1236 1136 void Tesselation::GoToFirst() const 1237 1137 { 1238 Info FunctionInfo(__func__);1138 Info FunctionInfo(__func__); 1239 1139 InternalPointer = PointsOnBoundary.begin(); 1240 } 1241 ; 1140 }; 1242 1141 1243 1142 /** PointCloud implementation of GoToLast. … … 1246 1145 void Tesselation::GoToLast() const 1247 1146 { 1248 Info FunctionInfo(__func__);1147 Info FunctionInfo(__func__); 1249 1148 InternalPointer = PointsOnBoundary.end(); 1250 1149 InternalPointer--; 1251 } 1252 ; 1150 }; 1253 1151 1254 1152 /** PointCloud implementation of IsEmpty. 1255 1153 * Uses PointsOnBoundary and STL stuff. 1256 */ 1154 */ 1257 1155 bool Tesselation::IsEmpty() const 1258 1156 { 1259 Info FunctionInfo(__func__);1157 Info FunctionInfo(__func__); 1260 1158 return (PointsOnBoundary.empty()); 1261 } 1262 ; 1159 }; 1263 1160 1264 1161 /** PointCloud implementation of IsLast. 1265 1162 * Uses PointsOnBoundary and STL stuff. 1266 */ 1163 */ 1267 1164 bool Tesselation::IsEnd() const 1268 1165 { 1269 Info FunctionInfo(__func__);1166 Info FunctionInfo(__func__); 1270 1167 return (InternalPointer == PointsOnBoundary.end()); 1271 } 1272 ; 1168 }; 1169 1273 1170 1274 1171 /** Gueses first starting triangle of the convex envelope. … … 1279 1176 void Tesselation::GuessStartingTriangle() 1280 1177 { 1281 Info FunctionInfo(__func__);1178 Info FunctionInfo(__func__); 1282 1179 // 4b. create a starting triangle 1283 1180 // 4b1. create all distances … … 1289 1186 1290 1187 // with A chosen, take each pair B,C and sort 1291 if (A != PointsOnBoundary.end()) { 1292 B = A; 1293 B++; 1294 for (; B != PointsOnBoundary.end(); B++) { 1295 C = B; 1296 C++; 1297 for (; C != PointsOnBoundary.end(); C++) { 1298 tmp = A->second->node->node->DistanceSquared(B->second->node->node); 1299 distance = tmp * tmp; 1300 tmp = A->second->node->node->DistanceSquared(C->second->node->node); 1301 distance += tmp * tmp; 1302 tmp = B->second->node->node->DistanceSquared(C->second->node->node); 1303 distance += tmp * tmp; 1304 DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C))); 1305 } 1306 } 1307 } 1188 if (A != PointsOnBoundary.end()) 1189 { 1190 B = A; 1191 B++; 1192 for (; B != PointsOnBoundary.end(); B++) 1193 { 1194 C = B; 1195 C++; 1196 for (; C != PointsOnBoundary.end(); C++) 1197 { 1198 tmp = A->second->node->node->DistanceSquared(B->second->node->node); 1199 distance = tmp * tmp; 1200 tmp = A->second->node->node->DistanceSquared(C->second->node->node); 1201 distance += tmp * tmp; 1202 tmp = B->second->node->node->DistanceSquared(C->second->node->node); 1203 distance += tmp * tmp; 1204 DistanceMMap.insert(DistanceMultiMapPair(distance, pair<PointMap::iterator, PointMap::iterator> (B, C))); 1205 } 1206 } 1207 } 1308 1208 // // listing distances 1309 1209 // Log() << Verbose(1) << "Listing DistanceMMap:"; … … 1315 1215 // 1. we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate 1316 1216 DistanceMultiMap::iterator baseline = DistanceMMap.begin(); 1317 for (; baseline != DistanceMMap.end(); baseline++) { 1318 // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate 1319 // 2. next, we have to check whether all points reside on only one side of the triangle 1320 // 3. construct plane vector 1321 PlaneVector.MakeNormalVector(A->second->node->node, baseline->second.first->second->node->node, baseline->second.second->second->node->node); 1322 DoLog(2) && (Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl); 1323 // 4. loop over all points 1324 double sign = 0.; 1325 PointMap::iterator checker = PointsOnBoundary.begin(); 1326 for (; checker != PointsOnBoundary.end(); checker++) { 1327 // (neglecting A,B,C) 1328 if ((checker == A) || (checker == baseline->second.first) || (checker == baseline->second.second)) 1329 continue; 1330 // 4a. project onto plane vector 1331 TrialVector.CopyVector(checker->second->node->node); 1332 TrialVector.SubtractVector(A->second->node->node); 1333 distance = TrialVector.ScalarProduct(&PlaneVector); 1334 if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok 1335 continue; 1336 DoLog(2) && (Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl); 1337 tmp = distance / fabs(distance); 1338 // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle) 1339 if ((sign != 0) && (tmp != sign)) { 1340 // 4c. If so, break 4. loop and continue with next candidate in 1. loop 1341 DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name << "," << baseline->second.second->second->node->Name << " leaves " << checker->second->node->Name << " outside the convex hull." << endl); 1342 break; 1343 } else { // note the sign for later 1344 DoLog(2) && (Log() << Verbose(2) << "Current candidates: " << A->second->node->Name << "," << baseline->second.first->second->node->Name << "," << baseline->second.second->second->node->Name << " leave " << checker->second->node->Name << " inside the convex hull." << endl); 1345 sign = tmp; 1346 } 1347 // 4d. Check whether the point is inside the triangle (check distance to each node 1348 tmp = checker->second->node->node->DistanceSquared(A->second->node->node); 1349 int innerpoint = 0; 1350 if ((tmp < A->second->node->node->DistanceSquared(baseline->second.first->second->node->node)) && (tmp < A->second->node->node->DistanceSquared(baseline->second.second->second->node->node))) 1351 innerpoint++; 1352 tmp = checker->second->node->node->DistanceSquared(baseline->second.first->second->node->node); 1353 if ((tmp < baseline->second.first->second->node->node->DistanceSquared(A->second->node->node)) && (tmp < baseline->second.first->second->node->node->DistanceSquared(baseline->second.second->second->node->node))) 1354 innerpoint++; 1355 tmp = checker->second->node->node->DistanceSquared(baseline->second.second->second->node->node); 1356 if ((tmp < baseline->second.second->second->node->node->DistanceSquared(baseline->second.first->second->node->node)) && (tmp < baseline->second.second->second->node->node->DistanceSquared(A->second->node->node))) 1357 innerpoint++; 1358 // 4e. If so, break 4. loop and continue with next candidate in 1. loop 1359 if (innerpoint == 3) 1360 break; 1361 } 1362 // 5. come this far, all on same side? Then break 1. loop and construct triangle 1363 if (checker == PointsOnBoundary.end()) { 1364 DoLog(2) && (Log() << Verbose(2) << "Looks like we have a candidate!" << endl); 1365 break; 1366 } 1367 } 1368 if (baseline != DistanceMMap.end()) { 1369 BPS[0] = baseline->second.first->second; 1370 BPS[1] = baseline->second.second->second; 1371 BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1372 BPS[0] = A->second; 1373 BPS[1] = baseline->second.second->second; 1374 BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1375 BPS[0] = baseline->second.first->second; 1376 BPS[1] = A->second; 1377 BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1378 1379 // 4b3. insert created triangle 1380 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 1381 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1382 TrianglesOnBoundaryCount++; 1383 for (int i = 0; i < NDIM; i++) { 1384 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i])); 1385 LinesOnBoundaryCount++; 1386 } 1387 1388 DoLog(1) && (Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl); 1389 } else { 1390 DoeLog(0) && (eLog() << Verbose(0) << "No starting triangle found." << endl); 1391 } 1217 for (; baseline != DistanceMMap.end(); baseline++) 1218 { 1219 // we take from the smallest sum of squared distance as the base line BC (with peak A) onward as the triangle candidate 1220 // 2. next, we have to check whether all points reside on only one side of the triangle 1221 // 3. construct plane vector 1222 PlaneVector.MakeNormalVector(A->second->node->node, 1223 baseline->second.first->second->node->node, 1224 baseline->second.second->second->node->node); 1225 Log() << Verbose(2) << "Plane vector of candidate triangle is " << PlaneVector << endl; 1226 // 4. loop over all points 1227 double sign = 0.; 1228 PointMap::iterator checker = PointsOnBoundary.begin(); 1229 for (; checker != PointsOnBoundary.end(); checker++) 1230 { 1231 // (neglecting A,B,C) 1232 if ((checker == A) || (checker == baseline->second.first) || (checker 1233 == baseline->second.second)) 1234 continue; 1235 // 4a. project onto plane vector 1236 TrialVector.CopyVector(checker->second->node->node); 1237 TrialVector.SubtractVector(A->second->node->node); 1238 distance = TrialVector.ScalarProduct(&PlaneVector); 1239 if (fabs(distance) < 1e-4) // we need to have a small epsilon around 0 which is still ok 1240 continue; 1241 Log() << Verbose(2) << "Projection of " << checker->second->node->Name << " yields distance of " << distance << "." << endl; 1242 tmp = distance / fabs(distance); 1243 // 4b. Any have different sign to than before? (i.e. would lie outside convex hull with this starting triangle) 1244 if ((sign != 0) && (tmp != sign)) 1245 { 1246 // 4c. If so, break 4. loop and continue with next candidate in 1. loop 1247 Log() << Verbose(2) << "Current candidates: " 1248 << A->second->node->Name << "," 1249 << baseline->second.first->second->node->Name << "," 1250 << baseline->second.second->second->node->Name << " leaves " 1251 << checker->second->node->Name << " outside the convex hull." 1252 << endl; 1253 break; 1254 } 1255 else 1256 { // note the sign for later 1257 Log() << Verbose(2) << "Current candidates: " 1258 << A->second->node->Name << "," 1259 << baseline->second.first->second->node->Name << "," 1260 << baseline->second.second->second->node->Name << " leave " 1261 << checker->second->node->Name << " inside the convex hull." 1262 << endl; 1263 sign = tmp; 1264 } 1265 // 4d. Check whether the point is inside the triangle (check distance to each node 1266 tmp = checker->second->node->node->DistanceSquared(A->second->node->node); 1267 int innerpoint = 0; 1268 if ((tmp < A->second->node->node->DistanceSquared( 1269 baseline->second.first->second->node->node)) && (tmp 1270 < A->second->node->node->DistanceSquared( 1271 baseline->second.second->second->node->node))) 1272 innerpoint++; 1273 tmp = checker->second->node->node->DistanceSquared( 1274 baseline->second.first->second->node->node); 1275 if ((tmp < baseline->second.first->second->node->node->DistanceSquared( 1276 A->second->node->node)) && (tmp 1277 < baseline->second.first->second->node->node->DistanceSquared( 1278 baseline->second.second->second->node->node))) 1279 innerpoint++; 1280 tmp = checker->second->node->node->DistanceSquared( 1281 baseline->second.second->second->node->node); 1282 if ((tmp < baseline->second.second->second->node->node->DistanceSquared( 1283 baseline->second.first->second->node->node)) && (tmp 1284 < baseline->second.second->second->node->node->DistanceSquared( 1285 A->second->node->node))) 1286 innerpoint++; 1287 // 4e. If so, break 4. loop and continue with next candidate in 1. loop 1288 if (innerpoint == 3) 1289 break; 1290 } 1291 // 5. come this far, all on same side? Then break 1. loop and construct triangle 1292 if (checker == PointsOnBoundary.end()) 1293 { 1294 Log() << Verbose(2) << "Looks like we have a candidate!" << endl; 1295 break; 1296 } 1297 } 1298 if (baseline != DistanceMMap.end()) 1299 { 1300 BPS[0] = baseline->second.first->second; 1301 BPS[1] = baseline->second.second->second; 1302 BLS[0] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1303 BPS[0] = A->second; 1304 BPS[1] = baseline->second.second->second; 1305 BLS[1] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1306 BPS[0] = baseline->second.first->second; 1307 BPS[1] = A->second; 1308 BLS[2] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1309 1310 // 4b3. insert created triangle 1311 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 1312 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1313 TrianglesOnBoundaryCount++; 1314 for (int i = 0; i < NDIM; i++) 1315 { 1316 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BTS->lines[i])); 1317 LinesOnBoundaryCount++; 1318 } 1319 1320 Log() << Verbose(1) << "Starting triangle is " << *BTS << "." << endl; 1321 } 1322 else 1323 { 1324 eLog() << Verbose(0) << "No starting triangle found." << endl; 1325 } 1392 1326 } 1393 1327 ; … … 1408 1342 void Tesselation::TesselateOnBoundary(const PointCloud * const cloud) 1409 1343 { 1410 Info FunctionInfo(__func__);1344 Info FunctionInfo(__func__); 1411 1345 bool flag; 1412 1346 PointMap::iterator winner; … … 1427 1361 // get peak point with respect to this base line's only triangle 1428 1362 BTS = baseline->second->triangles.begin()->second; // there is only one triangle so far 1429 DoLog(0) && (Log() << Verbose(0) << "Current baseline is between " << *(baseline->second) << "." << endl);1363 Log() << Verbose(0) << "Current baseline is between " << *(baseline->second) << "." << endl; 1430 1364 for (int i = 0; i < 3; i++) 1431 1365 if ((BTS->endpoints[i] != baseline->second->endpoints[0]) && (BTS->endpoints[i] != baseline->second->endpoints[1])) 1432 1366 peak = BTS->endpoints[i]; 1433 DoLog(1) && (Log() << Verbose(1) << " and has peak " << *peak << "." << endl);1367 Log() << Verbose(1) << " and has peak " << *peak << "." << endl; 1434 1368 1435 1369 // prepare some auxiliary vectors … … 1446 1380 CenterVector.AddVector(BTS->endpoints[i]->node->node); 1447 1381 CenterVector.Scale(1. / 3.); 1448 DoLog(2) && (Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl);1382 Log() << Verbose(2) << "CenterVector of base triangle is " << CenterVector << endl; 1449 1383 1450 1384 // normal vector of triangle … … 1453 1387 BTS->GetNormalVector(NormalVector); 1454 1388 NormalVector.CopyVector(&BTS->NormalVector); 1455 DoLog(2) && (Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl);1389 Log() << Verbose(2) << "NormalVector of base triangle is " << NormalVector << endl; 1456 1390 1457 1391 // vector in propagation direction (out of triangle) … … 1463 1397 if (PropagationVector.ScalarProduct(&TempVector) > 0) // make sure normal propagation vector points outward from baseline 1464 1398 PropagationVector.Scale(-1.); 1465 DoLog(2) && (Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl);1399 Log() << Verbose(2) << "PropagationVector of base triangle is " << PropagationVector << endl; 1466 1400 winner = PointsOnBoundary.end(); 1467 1401 … … 1469 1403 for (PointMap::iterator target = PointsOnBoundary.begin(); target != PointsOnBoundary.end(); target++) { 1470 1404 if ((target->second != baseline->second->endpoints[0]) && (target->second != baseline->second->endpoints[1])) { // don't take the same endpoints 1471 DoLog(1) && (Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl);1405 Log() << Verbose(1) << "Target point is " << *(target->second) << ":" << endl; 1472 1406 1473 1407 // first check direction, so that triangles don't intersect … … 1476 1410 VirtualNormalVector.ProjectOntoPlane(&NormalVector); 1477 1411 TempAngle = VirtualNormalVector.Angle(&PropagationVector); 1478 DoLog(2) && (Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl);1479 if (TempAngle > (M_PI /2.)) { // no bends bigger than Pi/2 (90 degrees)1480 DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl);1412 Log() << Verbose(2) << "VirtualNormalVector is " << VirtualNormalVector << " and PropagationVector is " << PropagationVector << "." << endl; 1413 if (TempAngle > (M_PI/2.)) { // no bends bigger than Pi/2 (90 degrees) 1414 Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", bad direction!" << endl; 1481 1415 continue; 1482 1416 } else 1483 DoLog(2) && (Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl);1417 Log() << Verbose(2) << "Angle on triangle plane between propagation direction and base line to " << *(target->second) << " is " << TempAngle << ", good direction!" << endl; 1484 1418 1485 1419 // check first and second endpoint (if any connecting line goes to target has at least not more than 1 triangle) … … 1487 1421 LineChecker[1] = baseline->second->endpoints[1]->lines.find(target->first); 1488 1422 if (((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[0]->second->triangles.size() == 2))) { 1489 DoLog(2) && (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);1423 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; 1490 1424 continue; 1491 1425 } 1492 1426 if (((LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (LineChecker[1]->second->triangles.size() == 2))) { 1493 DoLog(2) && (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);1427 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; 1494 1428 continue; 1495 1429 } … … 1497 1431 // check whether the envisaged triangle does not already exist (if both lines exist and have same endpoint) 1498 1432 if ((((LineChecker[0] != baseline->second->endpoints[0]->lines.end()) && (LineChecker[1] != baseline->second->endpoints[1]->lines.end()) && (GetCommonEndpoint(LineChecker[0]->second, LineChecker[1]->second) == peak)))) { 1499 DoLog(4) && (Log() << Verbose(4) << "Current target is peak!" << endl);1433 Log() << Verbose(4) << "Current target is peak!" << endl; 1500 1434 continue; 1501 1435 } … … 1508 1442 helper.ProjectOntoPlane(&TempVector); 1509 1443 if (fabs(helper.NormSquared()) < MYEPSILON) { 1510 DoLog(2) && (Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl);1444 Log() << Verbose(2) << "Chosen set of vectors is linear dependent." << endl; 1511 1445 continue; 1512 1446 } … … 1518 1452 TempVector.AddVector(baseline->second->endpoints[1]->node->node); 1519 1453 TempVector.AddVector(target->second->node->node); 1520 TempVector.Scale(1. /3.);1454 TempVector.Scale(1./3.); 1521 1455 TempVector.SubtractVector(Center); 1522 1456 // make it always point outward … … 1525 1459 // calculate angle 1526 1460 TempAngle = NormalVector.Angle(&VirtualNormalVector); 1527 DoLog(2) && (Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl);1461 Log() << Verbose(2) << "NormalVector is " << VirtualNormalVector << " and the angle is " << TempAngle << "." << endl; 1528 1462 if ((SmallestAngle - TempAngle) > MYEPSILON) { // set to new possible winner 1529 1463 SmallestAngle = TempAngle; 1530 1464 winner = target; 1531 DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);1465 Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl; 1532 1466 } else if (fabs(SmallestAngle - TempAngle) < MYEPSILON) { // check the angle to propagation, both possible targets are in one plane! (their normals have same angle) 1533 1467 // hence, check the angles to some normal direction from our base line but in this common plane of both targets... … … 1547 1481 SmallestAngle = TempAngle; 1548 1482 winner = target; 1549 DoLog(2) && (Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl);1483 Log() << Verbose(2) << "New winner " << *winner->second->node << " due to smaller angle " << TempAngle << " to propagation direction." << endl; 1550 1484 } else 1551 DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl);1485 Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle to propagation direction." << endl; 1552 1486 } else 1553 DoLog(2) && (Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl);1487 Log() << Verbose(2) << "Keeping old winner " << *winner->second->node << " due to smaller angle between normal vectors." << endl; 1554 1488 } 1555 1489 } // end of loop over all boundary points … … 1557 1491 // 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 1558 1492 if (winner != PointsOnBoundary.end()) { 1559 DoLog(0) && (Log() << Verbose(0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl);1493 Log() << Verbose(0) << "Winning target point is " << *(winner->second) << " with angle " << SmallestAngle << "." << endl; 1560 1494 // create the lins of not yet present 1561 1495 BLS[0] = baseline->second; … … 1587 1521 TrianglesOnBoundaryCount++; 1588 1522 } else { 1589 DoeLog(2) && (eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl);1523 eLog() << Verbose(2) << "I could not determine a winner for this baseline " << *(baseline->second) << "." << endl; 1590 1524 } 1591 1525 1592 1526 // 5d. If the set of lines is not yet empty, go to 5. and continue 1593 1527 } else 1594 DoLog(0) && (Log() << Verbose(0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl);1528 Log() << Verbose(0) << "Baseline candidate " << *(baseline->second) << " has a triangle count of " << baseline->second->triangles.size() << "." << endl; 1595 1529 } while (flag); 1596 1530 1597 1531 // exit 1598 delete (Center); 1599 } 1600 ; 1532 delete(Center); 1533 }; 1601 1534 1602 1535 /** Inserts all points outside of the tesselated surface into it by adding new triangles. … … 1608 1541 bool Tesselation::InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC) 1609 1542 { 1610 Info FunctionInfo(__func__);1543 Info FunctionInfo(__func__); 1611 1544 Vector Intersection, Normal; 1612 1545 TesselPoint *Walker = NULL; … … 1618 1551 cloud->GoToFirst(); 1619 1552 BoundaryPoints = new LinkedCell(this, 5.); 1620 while (!cloud->IsEnd()) { // we only have to go once through all points, as boundary can become only bigger1553 while (!cloud->IsEnd()) { // we only have to go once through all points, as boundary can become only bigger 1621 1554 if (AddFlag) { 1622 delete (BoundaryPoints);1555 delete(BoundaryPoints); 1623 1556 BoundaryPoints = new LinkedCell(this, 5.); 1624 1557 AddFlag = false; 1625 1558 } 1626 1559 Walker = cloud->GetPoint(); 1627 DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Walker << "." << endl);1560 Log() << Verbose(0) << "Current point is " << *Walker << "." << endl; 1628 1561 // get the next triangle 1629 1562 triangles = FindClosestTrianglesToVector(Walker->node, BoundaryPoints); 1630 1563 BTS = triangles->front(); 1631 1564 if ((triangles == NULL) || (BTS->ContainsBoundaryPoint(Walker))) { 1632 DoLog(0) && (Log() << Verbose(0) << "No triangles found, probably a tesselation point itself." << endl);1565 Log() << Verbose(0) << "No triangles found, probably a tesselation point itself." << endl; 1633 1566 cloud->GoToNext(); 1634 1567 continue; 1635 1568 } else { 1636 1569 } 1637 DoLog(0) && (Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl);1570 Log() << Verbose(0) << "Closest triangle is " << *BTS << "." << endl; 1638 1571 // get the intersection point 1639 1572 if (BTS->GetIntersectionInsideTriangle(Center, Walker->node, &Intersection)) { 1640 DoLog(0) && (Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl);1573 Log() << Verbose(0) << "We have an intersection at " << Intersection << "." << endl; 1641 1574 // we have the intersection, check whether in- or outside of boundary 1642 1575 if ((Center->DistanceSquared(Walker->node) - Center->DistanceSquared(&Intersection)) < -MYEPSILON) { 1643 1576 // inside, next! 1644 DoLog(0) && (Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl);1577 Log() << Verbose(0) << *Walker << " is inside wrt triangle " << *BTS << "." << endl; 1645 1578 } else { 1646 1579 // outside! 1647 DoLog(0) && (Log() << Verbose(0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl);1580 Log() << Verbose(0) << *Walker << " is outside wrt triangle " << *BTS << "." << endl; 1648 1581 class BoundaryLineSet *OldLines[3], *NewLines[3]; 1649 1582 class BoundaryPointSet *OldPoints[3], *NewPoint; 1650 1583 // store the three old lines and old points 1651 for (int i = 0; i < 3;i++) {1584 for (int i=0;i<3;i++) { 1652 1585 OldLines[i] = BTS->lines[i]; 1653 1586 OldPoints[i] = BTS->endpoints[i]; … … 1655 1588 Normal.CopyVector(&BTS->NormalVector); 1656 1589 // add Walker to boundary points 1657 DoLog(0) && (Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl);1590 Log() << Verbose(0) << "Adding " << *Walker << " to BoundaryPoints." << endl; 1658 1591 AddFlag = true; 1659 if (AddBoundaryPoint(Walker, 0))1592 if (AddBoundaryPoint(Walker,0)) 1660 1593 NewPoint = BPS[0]; 1661 1594 else 1662 1595 continue; 1663 1596 // remove triangle 1664 DoLog(0) && (Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl);1597 Log() << Verbose(0) << "Erasing triangle " << *BTS << "." << endl; 1665 1598 TrianglesOnBoundary.erase(BTS->Nr); 1666 delete (BTS);1599 delete(BTS); 1667 1600 // create three new boundary lines 1668 for (int i = 0; i < 3;i++) {1601 for (int i=0;i<3;i++) { 1669 1602 BPS[0] = NewPoint; 1670 1603 BPS[1] = OldPoints[i]; 1671 1604 NewLines[i] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); 1672 DoLog(1) && (Log() << Verbose(1) << "Creating new line " << *NewLines[i] << "." << endl);1605 Log() << Verbose(1) << "Creating new line " << *NewLines[i] << "." << endl; 1673 1606 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, NewLines[i])); // no need for check for unique insertion as BPS[0] is definitely a new one 1674 1607 LinesOnBoundaryCount++; 1675 1608 } 1676 1609 // create three new triangle with new point 1677 for (int i = 0; i < 3;i++) { // find all baselines1610 for (int i=0;i<3;i++) { // find all baselines 1678 1611 BLS[0] = OldLines[i]; 1679 1612 int n = 1; 1680 for (int j = 0; j < 3;j++) {1613 for (int j=0;j<3;j++) { 1681 1614 if (NewLines[j]->IsConnectedTo(BLS[0])) { 1682 if (n >2) {1683 DoeLog(2) && (eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl);1615 if (n>2) { 1616 eLog() << Verbose(2) << BLS[0] << " connects to all of the new lines?!" << endl; 1684 1617 return false; 1685 1618 } else … … 1692 1625 BTS->GetNormalVector(Normal); 1693 1626 Normal.Scale(-1.); 1694 DoLog(0) && (Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl);1627 Log() << Verbose(0) << "Created new triangle " << *BTS << "." << endl; 1695 1628 TrianglesOnBoundary.insert(TrianglePair(TrianglesOnBoundaryCount, BTS)); 1696 1629 TrianglesOnBoundaryCount++; … … 1698 1631 } 1699 1632 } else { // something is wrong with FindClosestTriangleToPoint! 1700 DoeLog(1) && (eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl);1633 eLog() << Verbose(1) << "The closest triangle did not produce an intersection!" << endl; 1701 1634 return false; 1702 1635 } … … 1705 1638 1706 1639 // exit 1707 delete (Center);1640 delete(Center); 1708 1641 return true; 1709 } 1710 ; 1642 }; 1711 1643 1712 1644 /** Adds a point to the tesselation::PointsOnBoundary list. … … 1717 1649 bool Tesselation::AddBoundaryPoint(TesselPoint * Walker, const int n) 1718 1650 { 1719 Info FunctionInfo(__func__);1651 Info FunctionInfo(__func__); 1720 1652 PointTestPair InsertUnique; 1721 1653 BPS[n] = new class BoundaryPointSet(Walker); … … 1725 1657 return true; 1726 1658 } else { 1727 delete (BPS[n]);1659 delete(BPS[n]); 1728 1660 BPS[n] = InsertUnique.first->second; 1729 1661 return false; … … 1739 1671 void Tesselation::AddTesselationPoint(TesselPoint* Candidate, const int n) 1740 1672 { 1741 Info FunctionInfo(__func__);1673 Info FunctionInfo(__func__); 1742 1674 PointTestPair InsertUnique; 1743 1675 TPS[n] = new class BoundaryPointSet(Candidate); … … 1747 1679 } else { 1748 1680 delete TPS[n]; 1749 DoLog(0) && (Log() << Verbose(0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl);1681 Log() << Verbose(0) << "Node " << *((InsertUnique.first)->second->node) << " is already present in PointsOnBoundary." << endl; 1750 1682 TPS[n] = (InsertUnique.first)->second; 1751 1683 } … … 1760 1692 void Tesselation::SetTesselationPoint(TesselPoint* Candidate, const int n) const 1761 1693 { 1762 Info FunctionInfo(__func__);1694 Info FunctionInfo(__func__); 1763 1695 PointMap::const_iterator FindPoint = PointsOnBoundary.find(Candidate->nr); 1764 1696 if (FindPoint != PointsOnBoundary.end()) … … 1766 1698 else 1767 1699 TPS[n] = NULL; 1768 } 1769 ; 1700 }; 1770 1701 1771 1702 /** Function tries to add line from current Points in BPS to BoundaryLineSet. 1772 1703 * If successful it raises the line count and inserts the new line into the BLS, 1773 1704 * if unsuccessful, it writes the line which had been present into the BLS, deleting the new constructed one. 1774 * @param *OptCenter desired OptCenter if there are more than one candidate line1775 * @param *candidate third point of the triangle to be, for checking between multiple open line candidates1776 1705 * @param *a first endpoint 1777 1706 * @param *b second endpoint 1778 1707 * @param n index of Tesselation::BLS giving the line with both endpoints 1779 1708 */ 1780 void Tesselation::AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) 1781 { 1709 void Tesselation::AddTesselationLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) { 1782 1710 bool insertNewLine = true; 1711 1783 1712 LineMap::iterator FindLine = a->lines.find(b->node->nr); 1784 BoundaryLineSet *WinningLine = NULL;1785 1713 if (FindLine != a->lines.end()) { 1786 DoLog(1) && (Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl);1787 1788 pair<LineMap::iterator, LineMap::iterator> FindPair;1714 Log() << Verbose(1) << "INFO: There is at least one line between " << *a << " and " << *b << ": " << *(FindLine->second) << "." << endl; 1715 1716 pair<LineMap::iterator,LineMap::iterator> FindPair; 1789 1717 FindPair = a->lines.equal_range(b->node->nr); 1790 1718 1791 for (FindLine = FindPair.first; (FindLine != FindPair.second) && (insertNewLine); FindLine++) { 1792 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl); 1719 for (FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) { 1793 1720 // If there is a line with less than two attached triangles, we don't need a new line. 1794 if (FindLine->second->triangles.size() == 1) {1795 CandidateMap::iterator Finder = OpenLines.find(FindLine->second);1796 if (!Finder->second->pointlist.empty())1797 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl); 1798 else1799 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate." << endl);1800 // get open line1801 for (TesselPointList::const_iterator CandidateChecker = Finder->second->pointlist.begin(); CandidateChecker != Finder->second->pointlist.end(); ++CandidateChecker) { 1802 if ((*(CandidateChecker) == candidate->node) && (OptCenter == NULL || OptCenter->DistanceSquared(&Finder->second->OptCenter) < MYEPSILON )) { // stop searching if candidate matches1803 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Candidate " << *(*CandidateChecker) << " has the right center " << Finder->second->OptCenter << "." << endl);1804 insertNewLine = false;1805 WinningLine = FindLine->second;1806 break;1807 } else {1808 DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *(*CandidateChecker) << "'s center " << Finder->second->OptCenter << " does not match desired on " << *OptCenter << "." << endl);1809 }1721 if (FindLine->second->triangles.size() < 2) { 1722 insertNewLine = false; 1723 Log() << Verbose(0) << "Using existing line " << *FindLine->second << endl; 1724 1725 BPS[0] = FindLine->second->endpoints[0]; 1726 BPS[1] = FindLine->second->endpoints[1]; 1727 BLS[n] = FindLine->second; 1728 1729 // remove existing line from OpenLines 1730 CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]); 1731 if (CandidateLine != OpenLines.end()) { 1732 Log() << Verbose(1) << " Removing line from OpenLines." << endl; 1733 delete(CandidateLine->second); 1734 OpenLines.erase(CandidateLine); 1735 } else { 1736 eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl; 1810 1737 } 1738 1739 break; 1811 1740 } 1812 1741 } … … 1814 1743 1815 1744 if (insertNewLine) { 1816 AddNewTesselationTriangleLine(a, b, n); 1817 } else { 1818 AddExistingTesselationTriangleLine(WinningLine, n); 1745 AlwaysAddTesselationTriangleLine(a, b, n); 1819 1746 } 1820 1747 } … … 1829 1756 * @param n index of Tesselation::BLS giving the line with both endpoints 1830 1757 */ 1831 void Tesselation::A ddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n)1832 { 1833 Info FunctionInfo(__func__);1834 DoLog(0) && (Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl);1758 void Tesselation::AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n) 1759 { 1760 Info FunctionInfo(__func__); 1761 Log() << Verbose(0) << "Adding open line [" << LinesOnBoundaryCount << "|" << *(a->node) << " and " << *(b->node) << "." << endl; 1835 1762 BPS[0] = a; 1836 1763 BPS[1] = b; 1837 BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); // this also adds the line to the local maps1764 BLS[n] = new class BoundaryLineSet(BPS, LinesOnBoundaryCount); // this also adds the line to the local maps 1838 1765 // add line to global map 1839 1766 LinesOnBoundary.insert(LinePair(LinesOnBoundaryCount, BLS[n])); … … 1842 1769 // also add to open lines 1843 1770 CandidateForTesselation *CFT = new CandidateForTesselation(BLS[n]); 1844 OpenLines.insert(pair<BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT)); 1845 } 1846 ; 1847 1848 /** Uses an existing line for a new triangle. 1849 * Sets Tesselation::BLS[\a n] and removes the lines from Tesselation::OpenLines. 1850 * \param *FindLine the line to add 1851 * \param n index of the line to set in Tesselation::BLS 1852 */ 1853 void Tesselation::AddExistingTesselationTriangleLine(class BoundaryLineSet *Line, int n) 1854 { 1855 Info FunctionInfo(__func__); 1856 DoLog(0) && (Log() << Verbose(0) << "Using existing line " << *Line << endl); 1857 1858 // set endpoints and line 1859 BPS[0] = Line->endpoints[0]; 1860 BPS[1] = Line->endpoints[1]; 1861 BLS[n] = Line; 1862 // remove existing line from OpenLines 1863 CandidateMap::iterator CandidateLine = OpenLines.find(BLS[n]); 1864 if (CandidateLine != OpenLines.end()) { 1865 DoLog(1) && (Log() << Verbose(1) << " Removing line from OpenLines." << endl); 1866 delete (CandidateLine->second); 1867 OpenLines.erase(CandidateLine); 1868 } else { 1869 DoeLog(1) && (eLog() << Verbose(1) << "Line exists and is attached to less than two triangles, but not in OpenLines!" << endl); 1870 } 1871 } 1872 ; 1771 OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (BLS[n], CFT)); 1772 }; 1873 1773 1874 1774 /** Function adds triangle to global list. … … 1877 1777 void Tesselation::AddTesselationTriangle() 1878 1778 { 1879 Info FunctionInfo(__func__);1880 DoLog(1) && (Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl);1779 Info FunctionInfo(__func__); 1780 Log() << Verbose(1) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1881 1781 1882 1782 // add triangle to global map … … 1888 1788 1889 1789 // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet 1890 } 1891 ; 1790 }; 1892 1791 1893 1792 /** Function adds triangle to global list. … … 1897 1796 void Tesselation::AddTesselationTriangle(const int nr) 1898 1797 { 1899 Info FunctionInfo(__func__);1900 DoLog(0) && (Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl);1798 Info FunctionInfo(__func__); 1799 Log() << Verbose(0) << "Adding triangle to global TrianglesOnBoundary map." << endl; 1901 1800 1902 1801 // add triangle to global map … … 1907 1806 1908 1807 // NOTE: add triangle to local maps is done in constructor of BoundaryTriangleSet 1909 } 1910 ; 1808 }; 1911 1809 1912 1810 /** Removes a triangle from the tesselation. … … 1917 1815 void Tesselation::RemoveTesselationTriangle(class BoundaryTriangleSet *triangle) 1918 1816 { 1919 Info FunctionInfo(__func__);1817 Info FunctionInfo(__func__); 1920 1818 if (triangle == NULL) 1921 1819 return; 1922 1820 for (int i = 0; i < 3; i++) { 1923 1821 if (triangle->lines[i] != NULL) { 1924 DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl);1822 Log() << Verbose(0) << "Removing triangle Nr." << triangle->Nr << " in line " << *triangle->lines[i] << "." << endl; 1925 1823 triangle->lines[i]->triangles.erase(triangle->Nr); 1926 1824 if (triangle->lines[i]->triangles.empty()) { 1927 DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl);1928 RemoveTesselationLine(triangle->lines[i]);1825 Log() << Verbose(0) << *triangle->lines[i] << " is no more attached to any triangle, erasing." << endl; 1826 RemoveTesselationLine(triangle->lines[i]); 1929 1827 } else { 1930 DoLog(0) && (Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: ");1931 OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL));1932 for (TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++)1933 DoLog(0) && (Log() << Verbose(0) << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t");1934 DoLog(0) && (Log() << Verbose(0) << endl);1935 // for (int j=0;j<2;j++) {1936 // Log() << Verbose(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": ";1937 // for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++)1938 // Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t";1939 // Log() << Verbose(0) << endl;1940 // }1828 Log() << Verbose(0) << *triangle->lines[i] << " is still attached to another triangle: "; 1829 OpenLines.insert(pair< BoundaryLineSet *, CandidateForTesselation *> (triangle->lines[i], NULL)); 1830 for(TriangleMap::iterator TriangleRunner = triangle->lines[i]->triangles.begin(); TriangleRunner != triangle->lines[i]->triangles.end(); TriangleRunner++) 1831 Log() << Verbose(0) << "[" << (TriangleRunner->second)->Nr << "|" << *((TriangleRunner->second)->endpoints[0]) << ", " << *((TriangleRunner->second)->endpoints[1]) << ", " << *((TriangleRunner->second)->endpoints[2]) << "] \t"; 1832 Log() << Verbose(0) << endl; 1833 // for (int j=0;j<2;j++) { 1834 // Log() << Verbose(0) << "Lines of endpoint " << *(triangle->lines[i]->endpoints[j]) << ": "; 1835 // for(LineMap::iterator LineRunner = triangle->lines[i]->endpoints[j]->lines.begin(); LineRunner != triangle->lines[i]->endpoints[j]->lines.end(); LineRunner++) 1836 // Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t"; 1837 // Log() << Verbose(0) << endl; 1838 // } 1941 1839 } 1942 triangle->lines[i] = NULL; // free'd or not: disconnect1840 triangle->lines[i] = NULL; // free'd or not: disconnect 1943 1841 } else 1944 DoeLog(1) && (eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl);1842 eLog() << Verbose(1) << "This line " << i << " has already been free'd." << endl; 1945 1843 } 1946 1844 1947 1845 if (TrianglesOnBoundary.erase(triangle->Nr)) 1948 DoLog(0) && (Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl); 1949 delete (triangle); 1950 } 1951 ; 1846 Log() << Verbose(0) << "Removing triangle Nr. " << triangle->Nr << "." << endl; 1847 delete(triangle); 1848 }; 1952 1849 1953 1850 /** Removes a line from the tesselation. … … 1957 1854 void Tesselation::RemoveTesselationLine(class BoundaryLineSet *line) 1958 1855 { 1959 Info FunctionInfo(__func__);1856 Info FunctionInfo(__func__); 1960 1857 int Numbers[2]; 1961 1858 … … 1978 1875 for (LineMap::iterator Runner = erasor.first; Runner != erasor.second; Runner++) 1979 1876 if ((*Runner).second == line) { 1980 DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);1877 Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl; 1981 1878 line->endpoints[i]->lines.erase(Runner); 1982 1879 break; … … 1984 1881 } else { // there's just a single line left 1985 1882 if (line->endpoints[i]->lines.erase(line->Nr)) 1986 DoLog(0) && (Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl);1883 Log() << Verbose(0) << "Removing Line Nr. " << line->Nr << " in boundary point " << *line->endpoints[i] << "." << endl; 1987 1884 } 1988 1885 if (line->endpoints[i]->lines.empty()) { 1989 DoLog(0) && (Log() << Verbose(0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl);1886 Log() << Verbose(0) << *line->endpoints[i] << " has no more lines it's attached to, erasing." << endl; 1990 1887 RemoveTesselationPoint(line->endpoints[i]); 1991 1888 } else { 1992 DoLog(0) && (Log() << Verbose(0) << *line->endpoints[i] << " has still lines it's attached to: ");1993 for (LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++)1994 DoLog(0) && (Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t");1995 DoLog(0) && (Log() << Verbose(0) << endl);1889 Log() << Verbose(0) << *line->endpoints[i] << " has still lines it's attached to: "; 1890 for(LineMap::iterator LineRunner = line->endpoints[i]->lines.begin(); LineRunner != line->endpoints[i]->lines.end(); LineRunner++) 1891 Log() << Verbose(0) << "[" << *(LineRunner->second) << "] \t"; 1892 Log() << Verbose(0) << endl; 1996 1893 } 1997 line->endpoints[i] = NULL; // free'd or not: disconnect1894 line->endpoints[i] = NULL; // free'd or not: disconnect 1998 1895 } else 1999 DoeLog(1) && (eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl);1896 eLog() << Verbose(1) << "Endpoint " << i << " has already been free'd." << endl; 2000 1897 } 2001 1898 if (!line->triangles.empty()) 2002 DoeLog(2) && (eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl);1899 eLog() << Verbose(2) << "Memory Leak! I " << *line << " am still connected to some triangles." << endl; 2003 1900 2004 1901 if (LinesOnBoundary.erase(line->Nr)) 2005 DoLog(0) && (Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl); 2006 delete (line); 2007 } 2008 ; 1902 Log() << Verbose(0) << "Removing line Nr. " << line->Nr << "." << endl; 1903 delete(line); 1904 }; 2009 1905 2010 1906 /** Removes a point from the tesselation. … … 2015 1911 void Tesselation::RemoveTesselationPoint(class BoundaryPointSet *point) 2016 1912 { 2017 Info FunctionInfo(__func__);1913 Info FunctionInfo(__func__); 2018 1914 if (point == NULL) 2019 1915 return; 2020 1916 if (PointsOnBoundary.erase(point->Nr)) 2021 DoLog(0) && (Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl); 2022 delete (point); 2023 } 2024 ; 2025 2026 /** Checks validity of a given sphere of a candidate line. 2027 * \sa CandidateForTesselation::CheckValidity(), which is more evolved. 2028 * We check CandidateForTesselation::OtherOptCenter 2029 * \param &CandidateLine contains other degenerated candidates which we have to subtract as well 2030 * \param RADIUS radius of sphere 2031 * \param *LC LinkedCell structure with other atoms 2032 * \return true - candidate triangle is degenerated, false - candidate triangle is not degenerated 2033 */ 2034 bool Tesselation::CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) const 2035 { 2036 Info FunctionInfo(__func__); 2037 2038 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking whether sphere contains no others points ..." << endl); 2039 bool flag = true; 2040 2041 DoLog(1) && (Log() << Verbose(1) << "Check by: draw sphere {" << CandidateLine.OtherOptCenter.x[0] << " " << CandidateLine.OtherOptCenter.x[1] << " " << CandidateLine.OtherOptCenter.x[2] << "} radius " << RADIUS << " resolution 30" << endl); 2042 // get all points inside the sphere 2043 TesselPointList *ListofPoints = LC->GetPointsInsideSphere(RADIUS, &CandidateLine.OtherOptCenter); 2044 2045 DoLog(1) && (Log() << Verbose(1) << "The following atoms are inside sphere at " << CandidateLine.OtherOptCenter << ":" << endl); 2046 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) 2047 DoLog(1) && (Log() << Verbose(1) << " " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl); 2048 2049 // remove triangles's endpoints 2050 for (int i = 0; i < 2; i++) 2051 ListofPoints->remove(CandidateLine.BaseLine->endpoints[i]->node); 2052 2053 // remove other candidates 2054 for (TesselPointList::const_iterator Runner = CandidateLine.pointlist.begin(); Runner != CandidateLine.pointlist.end(); ++Runner) 2055 ListofPoints->remove(*Runner); 2056 2057 // check for other points 2058 if (!ListofPoints->empty()) { 2059 DoLog(1) && (Log() << Verbose(1) << "CheckDegeneracy: There are still " << ListofPoints->size() << " points inside the sphere." << endl); 2060 flag = false; 2061 DoLog(1) && (Log() << Verbose(1) << "External atoms inside of sphere at " << CandidateLine.OtherOptCenter << ":" << endl); 2062 for (TesselPointList::const_iterator Runner = ListofPoints->begin(); Runner != ListofPoints->end(); ++Runner) 2063 DoLog(1) && (Log() << Verbose(1) << " " << *(*Runner) << " with distance " << (*Runner)->node->Distance(&CandidateLine.OtherOptCenter) << "." << endl); 2064 } 2065 delete (ListofPoints); 2066 2067 return flag; 2068 } 2069 ; 1917 Log() << Verbose(0) << "Removing point Nr. " << point->Nr << "." << endl; 1918 delete(point); 1919 }; 2070 1920 2071 1921 /** Checks whether the triangle consisting of the three points is already present. … … 2080 1930 int Tesselation::CheckPresenceOfTriangle(TesselPoint *Candidates[3]) const 2081 1931 { 2082 Info FunctionInfo(__func__);1932 Info FunctionInfo(__func__); 2083 1933 int adjacentTriangleCount = 0; 2084 1934 class BoundaryPointSet *Points[3]; … … 2102 1952 for (; (FindLine != Points[i]->lines.end()) && (FindLine->first == Points[j]->node->nr); FindLine++) { 2103 1953 TriangleMap *triangles = &FindLine->second->triangles; 2104 DoLog(1) && (Log() << Verbose(1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl);1954 Log() << Verbose(1) << "Current line is " << FindLine->first << ": " << *(FindLine->second) << " with triangles " << triangles << "." << endl; 2105 1955 for (TriangleMap::const_iterator FindTriangle = triangles->begin(); FindTriangle != triangles->end(); FindTriangle++) { 2106 1956 if (FindTriangle->second->IsPresentTupel(Points)) { … … 2108 1958 } 2109 1959 } 2110 DoLog(1) && (Log() << Verbose(1) << "end." << endl);1960 Log() << Verbose(1) << "end." << endl; 2111 1961 } 2112 1962 // Only one of the triangle lines must be considered for the triangle count. … … 2118 1968 } 2119 1969 2120 DoLog(0) && (Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl);1970 Log() << Verbose(0) << "Found " << adjacentTriangleCount << " adjacent triangles for the point set." << endl; 2121 1971 return adjacentTriangleCount; 2122 } 2123 ; 1972 }; 2124 1973 2125 1974 /** Checks whether the triangle consisting of the three points is already present. … … 2133 1982 class BoundaryTriangleSet * Tesselation::GetPresentTriangle(TesselPoint *Candidates[3]) 2134 1983 { 2135 Info FunctionInfo(__func__);1984 Info FunctionInfo(__func__); 2136 1985 class BoundaryTriangleSet *triangle = NULL; 2137 1986 class BoundaryPointSet *Points[3]; … … 2171 2020 2172 2021 return triangle; 2173 } 2174 ; 2022 }; 2023 2175 2024 2176 2025 /** Finds the starting triangle for FindNonConvexBorder(). … … 2181 2030 * \param RADIUS radius of virtual rolling sphere 2182 2031 * \param *LC LinkedCell structure with neighbouring TesselPoint's 2183 * \return true - a starting triangle has been created, false - no valid triple of points found 2184 */ 2185 bool Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC) 2186 { 2187 Info FunctionInfo(__func__); 2032 */ 2033 void Tesselation::FindStartingTriangle(const double RADIUS, const LinkedCell *LC) 2034 { 2035 Info FunctionInfo(__func__); 2188 2036 int i = 0; 2189 2037 TesselPoint* MaxPoint[NDIM]; 2190 2038 TesselPoint* Temporary; 2191 2039 double maxCoordinate[NDIM]; 2192 BoundaryLineSet *BaseLine = NULL;2040 BoundaryLineSet BaseLine; 2193 2041 Vector helper; 2194 2042 Vector Chord; 2195 2043 Vector SearchDirection; 2196 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers2044 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 2197 2045 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 2198 2046 Vector SphereCenter; … … 2207 2055 2208 2056 // 1. searching topmost point with respect to each axis 2209 for (int i = 0; i < NDIM;i++) { // each axis2210 LC->n[i] = LC->N[i] -1; // current axis is topmost cell2211 for (LC->n[(i + 1) % NDIM] = 0; LC->n[(i + 1) % NDIM] < LC->N[(i + 1) % NDIM]; LC->n[(i + 1) %NDIM]++)2212 for (LC->n[(i + 2) % NDIM] = 0; LC->n[(i + 2) % NDIM] < LC->N[(i + 2) % NDIM]; LC->n[(i + 2) %NDIM]++) {2213 const Linked Cell::LinkedNodes *List = LC->GetCurrentCell();2057 for (int i=0;i<NDIM;i++) { // each axis 2058 LC->n[i] = LC->N[i]-1; // current axis is topmost cell 2059 for (LC->n[(i+1)%NDIM]=0;LC->n[(i+1)%NDIM]<LC->N[(i+1)%NDIM];LC->n[(i+1)%NDIM]++) 2060 for (LC->n[(i+2)%NDIM]=0;LC->n[(i+2)%NDIM]<LC->N[(i+2)%NDIM];LC->n[(i+2)%NDIM]++) { 2061 const LinkedNodes *List = LC->GetCurrentCell(); 2214 2062 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 2215 2063 if (List != NULL) { 2216 for (Linked Cell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end();Runner++) {2064 for (LinkedNodes::const_iterator Runner = List->begin();Runner != List->end();Runner++) { 2217 2065 if ((*Runner)->node->x[i] > maxCoordinate[i]) { 2218 DoLog(1) && (Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl);2066 Log() << Verbose(1) << "New maximal for axis " << i << " node is " << *(*Runner) << " at " << *(*Runner)->node << "." << endl; 2219 2067 maxCoordinate[i] = (*Runner)->node->x[i]; 2220 2068 MaxPoint[i] = (*Runner); … … 2222 2070 } 2223 2071 } else { 2224 DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);2072 eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl; 2225 2073 } 2226 2074 } 2227 2075 } 2228 2076 2229 DoLog(1) && (Log() << Verbose(1) << "Found maximum coordinates: ");2230 for (int i = 0; i < NDIM;i++)2231 DoLog(0) && (Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t");2232 DoLog(0) && (Log() << Verbose(0) << endl);2077 Log() << Verbose(1) << "Found maximum coordinates: "; 2078 for (int i=0;i<NDIM;i++) 2079 Log() << Verbose(0) << i << ": " << *MaxPoint[i] << "\t"; 2080 Log() << Verbose(0) << endl; 2233 2081 2234 2082 BTS = NULL; 2235 for (int k = 0; k < NDIM;k++) {2083 for (int k=0;k<NDIM;k++) { 2236 2084 NormalVector.Zero(); 2237 2085 NormalVector.x[k] = 1.; 2238 BaseLine = new BoundaryLineSet(); 2239 BaseLine->endpoints[0] = new BoundaryPointSet(MaxPoint[k]); 2240 DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl); 2086 BaseLine.endpoints[0] = new BoundaryPointSet(MaxPoint[k]); 2087 Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl; 2241 2088 2242 2089 double ShortestAngle; 2243 2090 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. 2244 2091 2245 Temporary = NULL; 2246 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_... 2247 if (Temporary == NULL) { 2248 // have we found a second point? 2249 delete BaseLine; 2092 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_... 2093 if (Temporary == NULL) // have we found a second point? 2250 2094 continue; 2251 } 2252 BaseLine->endpoints[1] = new BoundaryPointSet(Temporary); 2095 BaseLine.endpoints[1] = new BoundaryPointSet(Temporary); 2253 2096 2254 2097 // construct center of circle 2255 CircleCenter.CopyVector(BaseLine ->endpoints[0]->node->node);2256 CircleCenter.AddVector(BaseLine ->endpoints[1]->node->node);2098 CircleCenter.CopyVector(BaseLine.endpoints[0]->node->node); 2099 CircleCenter.AddVector(BaseLine.endpoints[1]->node->node); 2257 2100 CircleCenter.Scale(0.5); 2258 2101 2259 2102 // construct normal vector of circle 2260 CirclePlaneNormal.CopyVector(BaseLine ->endpoints[0]->node->node);2261 CirclePlaneNormal.SubtractVector(BaseLine ->endpoints[1]->node->node);2103 CirclePlaneNormal.CopyVector(BaseLine.endpoints[0]->node->node); 2104 CirclePlaneNormal.SubtractVector(BaseLine.endpoints[1]->node->node); 2262 2105 2263 2106 double radius = CirclePlaneNormal.NormSquared(); 2264 double CircleRadius = sqrt(RADIUS * RADIUS - radius /4.);2107 double CircleRadius = sqrt(RADIUS*RADIUS - radius/4.); 2265 2108 2266 2109 NormalVector.ProjectOntoPlane(&CirclePlaneNormal); 2267 2110 NormalVector.Normalize(); 2268 ShortestAngle = 2. *M_PI; // This will indicate the quadrant.2111 ShortestAngle = 2.*M_PI; // This will indicate the quadrant. 2269 2112 2270 2113 SphereCenter.CopyVector(&NormalVector); … … 2274 2117 2275 2118 // look in one direction of baseline for initial candidate 2276 SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector); // whether we look "left" first or "right" first is not important ...2119 SearchDirection.MakeNormalVector(&CirclePlaneNormal, &NormalVector); // whether we look "left" first or "right" first is not important ... 2277 2120 2278 2121 // adding point 1 and point 2 and add the line between them 2279 DoLog(0) && (Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine->endpoints[0]->node << "." << endl);2280 DoLog(0) && (Log() << Verbose(0) << "Found second point is at " << *BaseLine->endpoints[1]->node << ".\n");2122 Log() << Verbose(0) << "Coordinates of start node at " << *BaseLine.endpoints[0]->node << "." << endl; 2123 Log() << Verbose(0) << "Found second point is at " << *BaseLine.endpoints[1]->node << ".\n"; 2281 2124 2282 2125 //Log() << Verbose(1) << "INFO: OldSphereCenter is at " << helper << ".\n"; 2283 CandidateForTesselation OptCandidates( BaseLine);2126 CandidateForTesselation OptCandidates(&BaseLine); 2284 2127 FindThirdPointForTesselation(NormalVector, SearchDirection, SphereCenter, OptCandidates, NULL, RADIUS, LC); 2285 DoLog(0) && (Log() << Verbose(0) << "List of third Points is:" << endl);2128 Log() << Verbose(0) << "List of third Points is:" << endl; 2286 2129 for (TesselPointList::iterator it = OptCandidates.pointlist.begin(); it != OptCandidates.pointlist.end(); it++) { 2287 DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl); 2288 } 2289 if (!OptCandidates.pointlist.empty()) { 2290 BTS = NULL; 2291 AddCandidatePolygon(OptCandidates, RADIUS, LC); 2292 } else { 2293 delete BaseLine; 2294 continue; 2295 } 2296 2297 if (BTS != NULL) { // we have created one starting triangle 2298 delete BaseLine; 2130 Log() << Verbose(0) << " " << *(*it) << endl; 2131 } 2132 2133 BTS = NULL; 2134 AddCandidateTriangle(OptCandidates); 2135 // delete(BaseLine.endpoints[0]); 2136 // delete(BaseLine.endpoints[1]); 2137 2138 if (BTS != NULL) // we have created one starting triangle 2299 2139 break; 2300 }else {2140 else { 2301 2141 // remove all candidates from the list and then the list itself 2302 2142 OptCandidates.pointlist.clear(); 2303 2143 } 2304 delete BaseLine; 2305 } 2306 2307 return (BTS != NULL); 2308 } 2309 ; 2144 } 2145 }; 2310 2146 2311 2147 /** Checks for a given baseline and a third point candidate whether baselines of the found triangle don't have even better candidates. … … 2378 2214 // if (fabs(OldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 2379 2215 // // rotated the wrong way! 2380 // DoeLog(1) && (eLog()<< Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);2216 // eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2381 2217 // } 2382 2218 // … … 2435 2271 // } 2436 2272 // } else { 2437 // DoeLog(2) && (eLog()<< Verbose(2) << "Baseline is connected to two triangles already?" << endl);2273 // eLog() << Verbose(2) << "Baseline is connected to two triangles already?" << endl; 2438 2274 // } 2439 2275 // } else { … … 2442 2278 // } 2443 2279 // } else { 2444 // DoeLog(1) && (eLog()<< Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl);2280 // eLog() << Verbose(1) << "Could not find the TesselPoint " << *ThirdNode << "." << endl; 2445 2281 // } 2446 2282 // … … 2456 2292 * @param *LC LinkedCell structure with neighbouring points 2457 2293 */ 2458 bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC) 2459 { 2460 Info FunctionInfo(__func__); 2294 bool Tesselation::FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC) 2295 { 2296 Info FunctionInfo(__func__); 2297 bool result = true; 2298 2461 2299 Vector CircleCenter; 2462 2300 Vector CirclePlaneNormal; … … 2464 2302 Vector SearchDirection; 2465 2303 Vector helper; 2466 BoundaryPointSet *ThirdPoint= NULL;2304 TesselPoint *ThirdNode = NULL; 2467 2305 LineMap::iterator testline; 2468 2306 double radius, CircleRadius; 2469 2307 2470 for (int i = 0; i < 3;i++)2471 if ((T.endpoints[i] != CandidateLine.BaseLine->endpoints[0]) && (T.endpoints[i] != CandidateLine.BaseLine->endpoints[1])) {2472 Third Point = T.endpoints[i];2308 for (int i=0;i<3;i++) 2309 if ((T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[0]->node) && (T.endpoints[i]->node != CandidateLine.BaseLine->endpoints[1]->node)) { 2310 ThirdNode = T.endpoints[i]->node; 2473 2311 break; 2474 2312 } 2475 DoLog(0) && (Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdPoint " << *ThirdPoint << " of triangle " << T << "." << endl); 2476 2477 CandidateLine.T = &T; 2313 Log() << Verbose(0) << "Current baseline is " << *CandidateLine.BaseLine << " with ThirdNode " << *ThirdNode << " of triangle " << T << "." << endl; 2478 2314 2479 2315 // construct center of circle … … 2488 2324 // calculate squared radius of circle 2489 2325 radius = CirclePlaneNormal.ScalarProduct(&CirclePlaneNormal); 2490 if (radius / 4. < RADIUS *RADIUS) {2326 if (radius/4. < RADIUS*RADIUS) { 2491 2327 // construct relative sphere center with now known CircleCenter 2492 2328 RelativeSphereCenter.CopyVector(&T.SphereCenter); 2493 2329 RelativeSphereCenter.SubtractVector(&CircleCenter); 2494 2330 2495 CircleRadius = RADIUS * RADIUS - radius /4.;2331 CircleRadius = RADIUS*RADIUS - radius/4.; 2496 2332 CirclePlaneNormal.Normalize(); 2497 DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);2498 2499 DoLog(1) && (Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl);2333 Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2334 2335 Log() << Verbose(1) << "INFO: OldSphereCenter is at " << T.SphereCenter << "." << endl; 2500 2336 2501 2337 // construct SearchDirection and an "outward pointer" 2502 2338 SearchDirection.MakeNormalVector(&RelativeSphereCenter, &CirclePlaneNormal); 2503 2339 helper.CopyVector(&CircleCenter); 2504 helper.SubtractVector(Third Point->node->node);2340 helper.SubtractVector(ThirdNode->node); 2505 2341 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)// ohoh, SearchDirection points inwards! 2506 2342 SearchDirection.Scale(-1.); 2507 DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);2343 Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2508 2344 if (fabs(RelativeSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { 2509 2345 // rotated the wrong way! 2510 DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl);2346 eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are still not orthogonal!" << endl; 2511 2347 } 2512 2348 2513 2349 // add third point 2514 FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, Third Point, RADIUS, LC);2350 FindThirdPointForTesselation(T.NormalVector, SearchDirection, T.SphereCenter, CandidateLine, ThirdNode, RADIUS, LC); 2515 2351 2516 2352 } else { 2517 DoLog(0) && (Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl);2353 Log() << Verbose(0) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and base triangle " << T << " is too big!" << endl; 2518 2354 } 2519 2355 2520 2356 if (CandidateLine.pointlist.empty()) { 2521 DoeLog(2) && (eLog() << Verbose(2) << "Could not find a suitable candidate." << endl);2357 eLog() << Verbose(2) << "Could not find a suitable candidate." << endl; 2522 2358 return false; 2523 2359 } 2524 DoLog(0) && (Log() << Verbose(0) << "Third Points are: " << endl);2360 Log() << Verbose(0) << "Third Points are: " << endl; 2525 2361 for (TesselPointList::iterator it = CandidateLine.pointlist.begin(); it != CandidateLine.pointlist.end(); ++it) { 2526 DoLog(0) && (Log() << Verbose(0) << " " << *(*it) << endl);2362 Log() << Verbose(0) << " " << *(*it) << endl; 2527 2363 } 2528 2364 2529 2365 return true; 2530 } 2531 ; 2532 2533 /** Walks through Tesselation::OpenLines() and finds candidates for newly created ones. 2534 * \param *&LCList atoms in LinkedCell list 2535 * \param RADIUS radius of the virtual sphere 2536 * \return true - for all open lines without candidates so far, a candidate has been found, 2537 * false - at least one open line without candidate still 2538 */ 2539 bool Tesselation::FindCandidatesforOpenLines(const double RADIUS, const LinkedCell *&LCList) 2540 { 2541 bool TesselationFailFlag = true; 2542 CandidateForTesselation *baseline = NULL; 2543 BoundaryTriangleSet *T = NULL; 2544 2545 for (CandidateMap::iterator Runner = OpenLines.begin(); Runner != OpenLines.end(); Runner++) { 2546 baseline = Runner->second; 2547 if (baseline->pointlist.empty()) { 2548 assert((baseline->BaseLine->triangles.size() == 1) && ("Open line without exactly one attached triangle")); 2549 T = (((baseline->BaseLine->triangles.begin()))->second); 2550 DoLog(1) && (Log() << Verbose(1) << "Finding best candidate for open line " << *baseline->BaseLine << " of triangle " << *T << endl); 2551 TesselationFailFlag = TesselationFailFlag && FindNextSuitableTriangle(*baseline, *T, RADIUS, LCList); //the line is there, so there is a triangle, but only one. 2552 } 2553 } 2554 return TesselationFailFlag; 2555 } 2556 ; 2366 2367 // BoundaryLineSet *BaseRay = CandidateLine.BaseLine; 2368 // for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2369 // Log() << Verbose(0) << "Third point candidate is " << *(*it)->point 2370 // << " with circumsphere's center at " << (*it)->OptCenter << "." << endl; 2371 // Log() << Verbose(0) << "Baseline is " << *BaseRay << endl; 2372 // 2373 // // check whether all edges of the new triangle still have space for one more triangle (i.e. TriangleCount <2) 2374 // TesselPoint *PointCandidates[3]; 2375 // PointCandidates[0] = (*it)->point; 2376 // PointCandidates[1] = BaseRay->endpoints[0]->node; 2377 // PointCandidates[2] = BaseRay->endpoints[1]->node; 2378 // int existentTrianglesCount = CheckPresenceOfTriangle(PointCandidates); 2379 // 2380 // BTS = NULL; 2381 // // check for present edges and whether we reach better candidates from them 2382 // //if (HasOtherBaselineBetterCandidate(BaseRay, (*it)->point, ShortestAngle, RADIUS, LC) ) { 2383 // if (0) { 2384 // result = false; 2385 // break; 2386 // } else { 2387 // // If there is no triangle, add it regularly. 2388 // if (existentTrianglesCount == 0) { 2389 // AddTesselationPoint((*it)->point, 0); 2390 // AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2391 // AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2392 // 2393 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const )TPS)) { 2394 // CandidateLine.point = (*it)->point; 2395 // CandidateLine.OptCenter.CopyVector(&((*it)->OptCenter)); 2396 // CandidateLine.OtherOptCenter.CopyVector(&((*it)->OtherOptCenter)); 2397 // CandidateLine.ShortestAngle = ShortestAngle; 2398 // } else { 2399 //// eLog() << Verbose(1) << "This triangle consisting of "; 2400 //// Log() << Verbose(0) << *(*it)->point << ", "; 2401 //// Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2402 //// Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2403 //// Log() << Verbose(0) << "exists and is not added, as it 0x80000000006fc150(does not seem helpful!" << endl; 2404 // result = false; 2405 // } 2406 // } else if ((existentTrianglesCount >= 1) && (existentTrianglesCount <= 3)) { // If there is a planar region within the structure, we need this triangle a second time. 2407 // AddTesselationPoint((*it)->point, 0); 2408 // AddTesselationPoint(BaseRay->endpoints[0]->node, 1); 2409 // AddTesselationPoint(BaseRay->endpoints[1]->node, 2); 2410 // 2411 // // 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) 2412 // // i.e. at least one of the three lines must be present with TriangleCount <= 1 2413 // if (CheckLineCriteriaForDegeneratedTriangle((const BoundaryPointSet ** const)TPS) || CandidateLine.BaseLine->skipped) { 2414 // CandidateLine.point = (*it)->point; 2415 // CandidateLine.OptCenter.CopyVector(&(*it)->OptCenter); 2416 // CandidateLine.OtherOptCenter.CopyVector(&(*it)->OtherOptCenter); 2417 // CandidateLine.ShortestAngle = ShortestAngle+2.*M_PI; 2418 // 2419 // } else { 2420 //// 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; 2421 // result = false; 2422 // } 2423 // } else { 2424 //// Log() << Verbose(1) << "This triangle consisting of "; 2425 //// Log() << Verbose(0) << *(*it)->point << ", "; 2426 //// Log() << Verbose(0) << *BaseRay->endpoints[0]->node << " and "; 2427 //// Log() << Verbose(0) << *BaseRay->endpoints[1]->node << " "; 2428 //// Log() << Verbose(0) << "is invalid!" << endl; 2429 // result = false; 2430 // } 2431 // } 2432 // 2433 // // set baseline to new ray from ref point (here endpoints[0]->node) to current candidate (here (*it)->point)) 2434 // BaseRay = BLS[0]; 2435 // if ((BTS != NULL) && (BTS->NormalVector.NormSquared() < MYEPSILON)) { 2436 // eLog() << Verbose(1) << "Triangle " << *BTS << " has zero normal vector!" << endl; 2437 // exit(255); 2438 // } 2439 // 2440 // } 2441 // 2442 // // remove all candidates from the list and then the list itself 2443 // class CandidateForTesselation *remover = NULL; 2444 // for (CandidateList::iterator it = OptCandidates->begin(); it != OptCandidates->end(); ++it) { 2445 // remover = *it; 2446 // delete(remover); 2447 // } 2448 // delete(OptCandidates); 2449 return result; 2450 }; 2557 2451 2558 2452 /** Adds the present line and candidate point from \a &CandidateLine to the Tesselation. 2559 2453 * \param CandidateLine triangle to add 2560 * \param RADIUS Radius of sphere 2561 * \param *LC LinkedCell structure 2562 * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in 2563 * AddTesselationLine() in AddCandidateTriangle() 2564 */ 2565 void Tesselation::AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC) 2566 { 2567 Info FunctionInfo(__func__); 2454 * \NOTE we need the copy operator here as the original CandidateForTesselation is removed in AddTesselationLine() 2455 */ 2456 void Tesselation::AddCandidateTriangle(CandidateForTesselation CandidateLine) 2457 { 2458 Info FunctionInfo(__func__); 2568 2459 Vector Center; 2569 2460 TesselPoint * const TurningPoint = CandidateLine.BaseLine->endpoints[0]->node; 2570 TesselPointList::iterator Runner;2571 TesselPointList::iterator Sprinter;2572 2461 2573 2462 // fill the set of neighbours … … 2578 2467 TesselPointList *connectedClosestPoints = GetCircleOfSetOfPoints(&SetOfNeighbours, TurningPoint, CandidateLine.BaseLine->endpoints[1]->node->node); 2579 2468 2580 DoLog(0) && (Log() << Verbose(0) << "List of Candidates for Turning Point " << *TurningPoint << ":" << endl); 2469 // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles) 2470 Log() << Verbose(0) << "List of Candidates for Turning Point: " << *TurningPoint << "." << endl; 2581 2471 for (TesselPointList::iterator TesselRunner = connectedClosestPoints->begin(); TesselRunner != connectedClosestPoints->end(); ++TesselRunner) 2582 DoLog(0) && (Log() << Verbose(0) << " " << **TesselRunner << endl); 2583 2584 // go through all angle-sorted candidates (in degenerate n-nodes case we may have to add multiple triangles) 2585 Runner = connectedClosestPoints->begin(); 2586 Sprinter = Runner; 2472 Log() << Verbose(0) << **TesselRunner << endl; 2473 TesselPointList::iterator Runner = connectedClosestPoints->begin(); 2474 TesselPointList::iterator Sprinter = Runner; 2587 2475 Sprinter++; 2588 while (Sprinter != connectedClosestPoints->end()) { 2589 DoLog(0) && (Log() << Verbose(0) << "Current Runner is " << *(*Runner) << " and sprinter is " << *(*Sprinter) << "." << endl); 2590 2476 while(Sprinter != connectedClosestPoints->end()) { 2477 // add the points 2591 2478 AddTesselationPoint(TurningPoint, 0); 2592 AddTesselationPoint(*Runner, 1); 2593 AddTesselationPoint(*Sprinter, 2); 2594 2595 AddCandidateTriangle(CandidateLine, Opt); 2596 2479 AddTesselationPoint((*Runner), 1); 2480 AddTesselationPoint((*Sprinter), 2); 2481 2482 // add the lines 2483 AddTesselationLine(TPS[0], TPS[1], 0); 2484 AddTesselationLine(TPS[0], TPS[2], 1); 2485 AddTesselationLine(TPS[1], TPS[2], 2); 2486 2487 // add the triangles 2488 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2489 AddTesselationTriangle(); 2490 BTS->GetCenter(&Center); 2491 Center.SubtractVector(&CandidateLine.OptCenter); 2492 BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter); 2493 BTS->GetNormalVector(Center); 2494 2495 Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << "." << endl; 2597 2496 Runner = Sprinter; 2598 2497 Sprinter++; 2599 if (Sprinter != connectedClosestPoints->end()) { 2600 // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked) 2601 FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OptCenter); // Assume BTS contains last triangle 2602 DoLog(0) && (Log() << Verbose(0) << " There are still more triangles to add." << endl); 2603 } 2604 // pick candidates for other open lines as well 2605 FindCandidatesforOpenLines(RADIUS, LC); 2606 2607 // check whether we add a degenerate or a normal triangle 2608 if (CheckDegeneracy(CandidateLine, RADIUS, LC)) { 2609 // add normal and degenerate triangles 2610 DoLog(1) && (Log() << Verbose(1) << "Triangle of endpoints " << *TPS[0] << "," << *TPS[1] << " and " << *TPS[2] << " is degenerated, adding both sides." << endl); 2611 AddCandidateTriangle(CandidateLine, OtherOpt); 2612 2613 if (Sprinter != connectedClosestPoints->end()) { 2614 // fill the internal open lines with its respective candidate (otherwise lines in degenerate case are not picked) 2615 FindDegeneratedCandidatesforOpenLines(*Sprinter, &CandidateLine.OtherOptCenter); 2616 } 2617 // pick candidates for other open lines as well 2618 FindCandidatesforOpenLines(RADIUS, LC); 2619 } 2620 } 2621 delete (connectedClosestPoints); 2622 }; 2623 2624 /** for polygons (multiple candidates for a baseline) sets internal edges to the correct next candidate. 2625 * \param *Sprinter next candidate to which internal open lines are set 2626 * \param *OptCenter OptCenter for this candidate 2627 */ 2628 void Tesselation::FindDegeneratedCandidatesforOpenLines(TesselPoint * const Sprinter, const Vector * const OptCenter) 2629 { 2630 Info FunctionInfo(__func__); 2631 2632 pair<LineMap::iterator, LineMap::iterator> FindPair = TPS[0]->lines.equal_range(TPS[2]->node->nr); 2633 for (LineMap::const_iterator FindLine = FindPair.first; FindLine != FindPair.second; FindLine++) { 2634 DoLog(1) && (Log() << Verbose(1) << "INFO: Checking line " << *(FindLine->second) << " ..." << endl); 2635 // If there is a line with less than two attached triangles, we don't need a new line. 2636 if (FindLine->second->triangles.size() == 1) { 2637 CandidateMap::iterator Finder = OpenLines.find(FindLine->second); 2638 if (!Finder->second->pointlist.empty()) 2639 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with candidate " << **(Finder->second->pointlist.begin()) << "." << endl); 2640 else { 2641 DoLog(1) && (Log() << Verbose(1) << "INFO: line " << *(FindLine->second) << " is open with no candidate, setting to next Sprinter" << (*Sprinter) << endl); 2642 Finder->second->T = BTS; // is last triangle 2643 Finder->second->pointlist.push_back(Sprinter); 2644 Finder->second->ShortestAngle = 0.; 2645 Finder->second->OptCenter.CopyVector(OptCenter); 2646 } 2647 } 2648 } 2649 }; 2650 2651 /** If a given \a *triangle is degenerated, this adds both sides. 2652 * i.e. the triangle with same BoundaryPointSet's but NormalVector in opposite direction. 2653 * Note that endpoints are stored in Tesselation::TPS 2654 * \param CandidateLine CanddiateForTesselation structure for the desired BoundaryLine 2655 * \param RADIUS radius of sphere 2656 * \param *LC pointer to LinkedCell structure 2657 */ 2658 void Tesselation::AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) 2659 { 2660 Info FunctionInfo(__func__); 2661 Vector Center; 2662 CandidateMap::const_iterator CandidateCheck = OpenLines.end(); 2663 BoundaryTriangleSet *triangle = NULL; 2664 2665 /// 1. Create or pick the lines for the first triangle 2666 DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for first triangle ..." << endl); 2667 for (int i = 0; i < 3; i++) { 2668 BLS[i] = NULL; 2669 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl); 2670 AddTesselationLine(&CandidateLine.OptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i); 2671 } 2672 2673 /// 2. create the first triangle and NormalVector and so on 2674 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding first triangle with center at " << CandidateLine.OptCenter << " ..." << endl); 2675 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2676 AddTesselationTriangle(); 2677 2678 // create normal vector 2679 BTS->GetCenter(&Center); 2680 Center.SubtractVector(&CandidateLine.OptCenter); 2681 BTS->SphereCenter.CopyVector(&CandidateLine.OptCenter); 2682 BTS->GetNormalVector(Center); 2683 // give some verbose output about the whole procedure 2684 if (CandidateLine.T != NULL) 2685 DoLog(0) && (Log() << Verbose(0) << "--> New triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl); 2686 else 2687 DoLog(0) && (Log() << Verbose(0) << "--> New starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl); 2688 triangle = BTS; 2689 2690 /// 3. Gather candidates for each new line 2691 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding candidates to new lines ..." << endl); 2692 for (int i = 0; i < 3; i++) { 2693 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl); 2694 CandidateCheck = OpenLines.find(BLS[i]); 2695 if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) { 2696 if (CandidateCheck->second->T == NULL) 2697 CandidateCheck->second->T = triangle; 2698 FindNextSuitableTriangle(*(CandidateCheck->second), *CandidateCheck->second->T, RADIUS, LC); 2699 } 2700 } 2701 2702 /// 4. Create or pick the lines for the second triangle 2703 DoLog(0) && (Log() << Verbose(0) << "INFO: Creating/Picking lines for second triangle ..." << endl); 2704 for (int i = 0; i < 3; i++) { 2705 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl); 2706 AddTesselationLine(&CandidateLine.OtherOptCenter, TPS[(i + 2) % 3], TPS[(i + 0) % 3], TPS[(i + 1) % 3], i); 2707 } 2708 2709 /// 5. create the second triangle and NormalVector and so on 2710 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangle with center at " << CandidateLine.OtherOptCenter << " ..." << endl); 2711 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2712 AddTesselationTriangle(); 2713 2714 BTS->SphereCenter.CopyVector(&CandidateLine.OtherOptCenter); 2715 // create normal vector in other direction 2716 BTS->GetNormalVector(&triangle->NormalVector); 2717 BTS->NormalVector.Scale(-1.); 2718 // give some verbose output about the whole procedure 2719 if (CandidateLine.T != NULL) 2720 DoLog(0) && (Log() << Verbose(0) << "--> New degenerate triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl); 2721 else 2722 DoLog(0) && (Log() << Verbose(0) << "--> New degenerate starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl); 2723 2724 /// 6. Adding triangle to new lines 2725 DoLog(0) && (Log() << Verbose(0) << "INFO: Adding second triangles to new lines ..." << endl); 2726 for (int i = 0; i < 3; i++) { 2727 DoLog(0) && (Log() << Verbose(0) << "Current line is between " << *TPS[(i + 0) % 3] << " and " << *TPS[(i + 1) % 3] << ":" << endl); 2728 CandidateCheck = OpenLines.find(BLS[i]); 2729 if ((CandidateCheck != OpenLines.end()) && (CandidateCheck->second->pointlist.empty())) { 2730 if (CandidateCheck->second->T == NULL) 2731 CandidateCheck->second->T = BTS; 2732 } 2733 } 2734 } 2735 ; 2736 2737 /** Adds a triangle to the Tesselation structure from three given TesselPoint's. 2738 * Note that endpoints are in Tesselation::TPS. 2739 * \param CandidateLine CandidateForTesselation structure contains other information 2740 * \param type which opt center to add (i.e. which side) and thus which NormalVector to take 2741 */ 2742 void Tesselation::AddCandidateTriangle(CandidateForTesselation &CandidateLine, enum centers type) 2743 { 2744 Info FunctionInfo(__func__); 2745 Vector Center; 2746 Vector *OptCenter = (type == Opt) ? &CandidateLine.OptCenter : &CandidateLine.OtherOptCenter; 2747 2748 // add the lines 2749 AddTesselationLine(OptCenter, TPS[2], TPS[0], TPS[1], 0); 2750 AddTesselationLine(OptCenter, TPS[1], TPS[0], TPS[2], 1); 2751 AddTesselationLine(OptCenter, TPS[0], TPS[1], TPS[2], 2); 2752 2753 // add the triangles 2754 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 2755 AddTesselationTriangle(); 2756 2757 // create normal vector 2758 BTS->GetCenter(&Center); 2759 Center.SubtractVector(OptCenter); 2760 BTS->SphereCenter.CopyVector(OptCenter); 2761 BTS->GetNormalVector(Center); 2762 2763 // give some verbose output about the whole procedure 2764 if (CandidateLine.T != NULL) 2765 DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "triangle with " << *BTS << " and normal vector " << BTS->NormalVector << ", from " << *CandidateLine.T << " and angle " << CandidateLine.ShortestAngle << "." << endl); 2766 else 2767 DoLog(0) && (Log() << Verbose(0) << "--> New" << ((type == OtherOpt) ? " degenerate " : " ") << "starting triangle with " << *BTS << " and normal vector " << BTS->NormalVector << " and no top triangle." << endl); 2768 } 2769 ; 2498 Log() << Verbose(0) << "Current Runner is " << **Runner << "." << endl; 2499 if (Sprinter != connectedClosestPoints->end()) 2500 Log() << Verbose(0) << " There are still more triangles to add." << endl; 2501 } 2502 delete(connectedClosestPoints); 2503 }; 2770 2504 2771 2505 /** Checks whether the quadragon of the two triangles connect to \a *Base is convex. … … 2778 2512 class BoundaryPointSet *Tesselation::IsConvexRectangle(class BoundaryLineSet *Base) 2779 2513 { 2780 Info FunctionInfo(__func__);2514 Info FunctionInfo(__func__); 2781 2515 class BoundaryPointSet *Spot = NULL; 2782 2516 class BoundaryLineSet *OtherBase; 2783 2517 Vector *ClosestPoint; 2784 2518 2785 int m =0;2786 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)2787 for (int j = 0; j < 3;j++) // all of their endpoints and baselines2519 int m=0; 2520 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2521 for (int j=0;j<3;j++) // all of their endpoints and baselines 2788 2522 if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints 2789 2523 BPS[m++] = runner->second->endpoints[j]; 2790 OtherBase = new class BoundaryLineSet(BPS, -1);2791 2792 DoLog(1) && (Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl);2793 DoLog(1) && (Log() << Verbose(1) << "INFO: Other base line is " << *OtherBase << "." << endl);2524 OtherBase = new class BoundaryLineSet(BPS,-1); 2525 2526 Log() << Verbose(1) << "INFO: Current base line is " << *Base << "." << endl; 2527 Log() << Verbose(1) << "INFO: Other base line is " << *OtherBase << "." << endl; 2794 2528 2795 2529 // get the closest point on each line to the other line … … 2797 2531 2798 2532 // delete the temporary other base line 2799 delete (OtherBase);2533 delete(OtherBase); 2800 2534 2801 2535 // get the distance vector from Base line to OtherBase line … … 2804 2538 BaseLine.CopyVector(Base->endpoints[1]->node->node); 2805 2539 BaseLine.SubtractVector(Base->endpoints[0]->node->node); 2806 for (int i = 0; i < 2;i++) {2540 for (int i=0;i<2;i++) { 2807 2541 DistanceToIntersection[i].CopyVector(ClosestPoint); 2808 2542 DistanceToIntersection[i].SubtractVector(Base->endpoints[i]->node->node); 2809 2543 distance[i] = BaseLine.ScalarProduct(&DistanceToIntersection[i]); 2810 2544 } 2811 delete (ClosestPoint);2812 if ((distance[0] * distance[1]) > 0) { // have same sign?2813 DoLog(1) && (Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl);2545 delete(ClosestPoint); 2546 if ((distance[0] * distance[1]) > 0) { // have same sign? 2547 Log() << Verbose(1) << "REJECT: Both SKPs have same sign: " << distance[0] << " and " << distance[1] << ". " << *Base << "' rectangle is concave." << endl; 2814 2548 if (distance[0] < distance[1]) { 2815 2549 Spot = Base->endpoints[0]; … … 2818 2552 } 2819 2553 return Spot; 2820 } else { // different sign, i.e. we are in between2821 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl);2554 } else { // different sign, i.e. we are in between 2555 Log() << Verbose(0) << "ACCEPT: Rectangle of triangles of base line " << *Base << " is convex." << endl; 2822 2556 return NULL; 2823 2557 } 2824 2558 2825 } 2826 ; 2559 }; 2827 2560 2828 2561 void Tesselation::PrintAllBoundaryPoints(ofstream *out) const 2829 2562 { 2830 Info FunctionInfo(__func__);2563 Info FunctionInfo(__func__); 2831 2564 // print all lines 2832 DoLog(0) && (Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl); 2833 for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin(); PointRunner != PointsOnBoundary.end(); PointRunner++) 2834 DoLog(0) && (Log() << Verbose(0) << *(PointRunner->second) << endl); 2835 } 2836 ; 2565 Log() << Verbose(0) << "Printing all boundary points for debugging:" << endl; 2566 for (PointMap::const_iterator PointRunner = PointsOnBoundary.begin();PointRunner != PointsOnBoundary.end(); PointRunner++) 2567 Log() << Verbose(0) << *(PointRunner->second) << endl; 2568 }; 2837 2569 2838 2570 void Tesselation::PrintAllBoundaryLines(ofstream *out) const 2839 2571 { 2840 Info FunctionInfo(__func__);2572 Info FunctionInfo(__func__); 2841 2573 // print all lines 2842 DoLog(0) && (Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl);2574 Log() << Verbose(0) << "Printing all boundary lines for debugging:" << endl; 2843 2575 for (LineMap::const_iterator LineRunner = LinesOnBoundary.begin(); LineRunner != LinesOnBoundary.end(); LineRunner++) 2844 DoLog(0) && (Log() << Verbose(0) << *(LineRunner->second) << endl); 2845 } 2846 ; 2576 Log() << Verbose(0) << *(LineRunner->second) << endl; 2577 }; 2847 2578 2848 2579 void Tesselation::PrintAllBoundaryTriangles(ofstream *out) const 2849 2580 { 2850 Info FunctionInfo(__func__);2581 Info FunctionInfo(__func__); 2851 2582 // print all triangles 2852 DoLog(0) && (Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl);2583 Log() << Verbose(0) << "Printing all boundary triangles for debugging:" << endl; 2853 2584 for (TriangleMap::const_iterator TriangleRunner = TrianglesOnBoundary.begin(); TriangleRunner != TrianglesOnBoundary.end(); TriangleRunner++) 2854 DoLog(0) && (Log() << Verbose(0) << *(TriangleRunner->second) << endl); 2855 } 2856 ; 2585 Log() << Verbose(0) << *(TriangleRunner->second) << endl; 2586 }; 2857 2587 2858 2588 /** For a given boundary line \a *Base and its two triangles, picks the central baseline that is "higher". … … 2863 2593 double Tesselation::PickFarthestofTwoBaselines(class BoundaryLineSet *Base) 2864 2594 { 2865 Info FunctionInfo(__func__);2595 Info FunctionInfo(__func__); 2866 2596 class BoundaryLineSet *OtherBase; 2867 2597 Vector *ClosestPoint[2]; 2868 2598 double volume; 2869 2599 2870 int m =0;2871 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)2872 for (int j = 0; j < 3;j++) // all of their endpoints and baselines2600 int m=0; 2601 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2602 for (int j=0;j<3;j++) // all of their endpoints and baselines 2873 2603 if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) // and neither of its endpoints 2874 2604 BPS[m++] = runner->second->endpoints[j]; 2875 OtherBase = new class BoundaryLineSet(BPS, -1);2876 2877 DoLog(0) && (Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl);2878 DoLog(0) && (Log() << Verbose(0) << "INFO: Other base line is " << *OtherBase << "." << endl);2605 OtherBase = new class BoundaryLineSet(BPS,-1); 2606 2607 Log() << Verbose(0) << "INFO: Current base line is " << *Base << "." << endl; 2608 Log() << Verbose(0) << "INFO: Other base line is " << *OtherBase << "." << endl; 2879 2609 2880 2610 // get the closest point on each line to the other line … … 2891 2621 2892 2622 // delete the temporary other base line and the closest points 2893 delete (ClosestPoint[0]);2894 delete (ClosestPoint[1]);2895 delete (OtherBase);2623 delete(ClosestPoint[0]); 2624 delete(ClosestPoint[1]); 2625 delete(OtherBase); 2896 2626 2897 2627 if (Distance.NormSquared() < MYEPSILON) { // check for intersection 2898 DoLog(0) && (Log() << Verbose(0) << "REJECT: Both lines have an intersection: Nothing to do." << endl);2628 Log() << Verbose(0) << "REJECT: Both lines have an intersection: Nothing to do." << endl; 2899 2629 return false; 2900 2630 } else { // check for sign against BaseLineNormal … … 2902 2632 BaseLineNormal.Zero(); 2903 2633 if (Base->triangles.size() < 2) { 2904 DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);2634 eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl; 2905 2635 return 0.; 2906 2636 } 2907 2637 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2908 DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);2638 Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl; 2909 2639 BaseLineNormal.AddVector(&(runner->second->NormalVector)); 2910 2640 } 2911 BaseLineNormal.Scale(1. /2.);2641 BaseLineNormal.Scale(1./2.); 2912 2642 2913 2643 if (Distance.ScalarProduct(&BaseLineNormal) > MYEPSILON) { // Distance points outwards, hence OtherBase higher than Base -> flip 2914 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl);2644 Log() << Verbose(0) << "ACCEPT: Other base line would be higher: Flipping baseline." << endl; 2915 2645 // calculate volume summand as a general tetraeder 2916 2646 return volume; 2917 } else { // Base higher than OtherBase -> do nothing2918 DoLog(0) && (Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl);2647 } else { // Base higher than OtherBase -> do nothing 2648 Log() << Verbose(0) << "REJECT: Base line is higher: Nothing to do." << endl; 2919 2649 return 0.; 2920 2650 } 2921 2651 } 2922 } 2923 ; 2652 }; 2924 2653 2925 2654 /** For a given baseline and its two connected triangles, flips the baseline. … … 2932 2661 class BoundaryLineSet * Tesselation::FlipBaseline(class BoundaryLineSet *Base) 2933 2662 { 2934 Info FunctionInfo(__func__);2663 Info FunctionInfo(__func__); 2935 2664 class BoundaryLineSet *OldLines[4], *NewLine; 2936 2665 class BoundaryPointSet *OldPoints[2]; 2937 2666 Vector BaseLineNormal; 2938 2667 int OldTriangleNrs[2], OldBaseLineNr; 2939 int i, m;2668 int i,m; 2940 2669 2941 2670 // calculate NormalVector for later use 2942 2671 BaseLineNormal.Zero(); 2943 2672 if (Base->triangles.size() < 2) { 2944 DoeLog(1) && (eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl);2673 eLog() << Verbose(1) << "Less than two triangles are attached to this baseline!" << endl; 2945 2674 return NULL; 2946 2675 } 2947 2676 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2948 DoLog(1) && (Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl);2677 Log() << Verbose(1) << "INFO: Adding NormalVector " << runner->second->NormalVector << " of triangle " << *(runner->second) << "." << endl; 2949 2678 BaseLineNormal.AddVector(&(runner->second->NormalVector)); 2950 2679 } 2951 BaseLineNormal.Scale(-1. /2.); // has to point inside for BoundaryTriangleSet::GetNormalVector()2680 BaseLineNormal.Scale(-1./2.); // has to point inside for BoundaryTriangleSet::GetNormalVector() 2952 2681 2953 2682 // get the two triangles 2954 2683 // gather four endpoints and four lines 2955 for (int j = 0; j < 4;j++)2684 for (int j=0;j<4;j++) 2956 2685 OldLines[j] = NULL; 2957 for (int j = 0; j < 2;j++)2686 for (int j=0;j<2;j++) 2958 2687 OldPoints[j] = NULL; 2959 i =0;2960 m =0;2961 DoLog(0) && (Log() << Verbose(0) << "The four old lines are: ");2962 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)2963 for (int j = 0; j < 3;j++) // all of their endpoints and baselines2688 i=0; 2689 m=0; 2690 Log() << Verbose(0) << "The four old lines are: "; 2691 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2692 for (int j=0;j<3;j++) // all of their endpoints and baselines 2964 2693 if (runner->second->lines[j] != Base) { // pick not the central baseline 2965 2694 OldLines[i++] = runner->second->lines[j]; 2966 DoLog(0) && (Log() << Verbose(0) << *runner->second->lines[j] << "\t");2695 Log() << Verbose(0) << *runner->second->lines[j] << "\t"; 2967 2696 } 2968 DoLog(0) && (Log() << Verbose(0) << endl);2969 DoLog(0) && (Log() << Verbose(0) << "The two old points are: ");2970 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++)2971 for (int j = 0; j < 3;j++) // all of their endpoints and baselines2697 Log() << Verbose(0) << endl; 2698 Log() << Verbose(0) << "The two old points are: "; 2699 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) 2700 for (int j=0;j<3;j++) // all of their endpoints and baselines 2972 2701 if (!Base->ContainsBoundaryPoint(runner->second->endpoints[j])) { // and neither of its endpoints 2973 2702 OldPoints[m++] = runner->second->endpoints[j]; 2974 DoLog(0) && (Log() << Verbose(0) << *runner->second->endpoints[j] << "\t");2703 Log() << Verbose(0) << *runner->second->endpoints[j] << "\t"; 2975 2704 } 2976 DoLog(0) && (Log() << Verbose(0) << endl);2705 Log() << Verbose(0) << endl; 2977 2706 2978 2707 // check whether everything is in place to create new lines and triangles 2979 if (i <4) {2980 DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);2708 if (i<4) { 2709 eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl; 2981 2710 return NULL; 2982 2711 } 2983 for (int j = 0; j < 4;j++)2712 for (int j=0;j<4;j++) 2984 2713 if (OldLines[j] == NULL) { 2985 DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl);2714 eLog() << Verbose(1) << "We have not gathered enough baselines!" << endl; 2986 2715 return NULL; 2987 2716 } 2988 for (int j = 0; j < 2;j++)2717 for (int j=0;j<2;j++) 2989 2718 if (OldPoints[j] == NULL) { 2990 DoeLog(1) && (eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl);2719 eLog() << Verbose(1) << "We have not gathered enough endpoints!" << endl; 2991 2720 return NULL; 2992 2721 } 2993 2722 2994 2723 // remove triangles and baseline removes itself 2995 DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting baseline " << *Base << " from global list." << endl);2724 Log() << Verbose(0) << "INFO: Deleting baseline " << *Base << " from global list." << endl; 2996 2725 OldBaseLineNr = Base->Nr; 2997 m =0;2998 for (TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) {2999 DoLog(0) && (Log() << Verbose(0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl);2726 m=0; 2727 for(TriangleMap::iterator runner = Base->triangles.begin(); runner != Base->triangles.end(); runner++) { 2728 Log() << Verbose(0) << "INFO: Deleting triangle " << *(runner->second) << "." << endl; 3000 2729 OldTriangleNrs[m++] = runner->second->Nr; 3001 2730 RemoveTesselationTriangle(runner->second); … … 3007 2736 NewLine = new class BoundaryLineSet(BPS, OldBaseLineNr); 3008 2737 LinesOnBoundary.insert(LinePair(OldBaseLineNr, NewLine)); // no need for check for unique insertion as NewLine is definitely a new one 3009 DoLog(0) && (Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl);2738 Log() << Verbose(0) << "INFO: Created new baseline " << *NewLine << "." << endl; 3010 2739 3011 2740 // construct new triangles with flipped baseline 3012 i =-1;2741 i=-1; 3013 2742 if (OldLines[0]->IsConnectedTo(OldLines[2])) 3014 i =2;2743 i=2; 3015 2744 if (OldLines[0]->IsConnectedTo(OldLines[3])) 3016 i =3;3017 if (i !=-1) {2745 i=3; 2746 if (i!=-1) { 3018 2747 BLS[0] = OldLines[0]; 3019 2748 BLS[1] = OldLines[i]; … … 3022 2751 BTS->GetNormalVector(BaseLineNormal); 3023 2752 AddTesselationTriangle(OldTriangleNrs[0]); 3024 DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);3025 3026 BLS[0] = (i ==2 ? OldLines[3] : OldLines[2]);2753 Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl; 2754 2755 BLS[0] = (i==2 ? OldLines[3] : OldLines[2]); 3027 2756 BLS[1] = OldLines[1]; 3028 2757 BLS[2] = NewLine; … … 3030 2759 BTS->GetNormalVector(BaseLineNormal); 3031 2760 AddTesselationTriangle(OldTriangleNrs[1]); 3032 DoLog(0) && (Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl);2761 Log() << Verbose(0) << "INFO: Created new triangle " << *BTS << "." << endl; 3033 2762 } else { 3034 DoeLog(0) && (eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl);2763 eLog() << Verbose(0) << "The four old lines do not connect, something's utterly wrong here!" << endl; 3035 2764 return NULL; 3036 2765 } 3037 2766 3038 2767 return NewLine; 3039 } 3040 ; 2768 }; 2769 3041 2770 3042 2771 /** Finds the second point of starting triangle. … … 3050 2779 void Tesselation::FindSecondPointForTesselation(TesselPoint* a, Vector Oben, TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC) 3051 2780 { 3052 Info FunctionInfo(__func__);2781 Info FunctionInfo(__func__); 3053 2782 Vector AngleCheck; 3054 2783 class TesselPoint* Candidate = NULL; … … 3059 2788 int Nupper[NDIM]; 3060 2789 3061 if (LC->SetIndexToNode(a)) { // get cell for the starting point3062 for (int i = 0; i < NDIM;i++) // store indices of this cell2790 if (LC->SetIndexToNode(a)) { // get cell for the starting point 2791 for(int i=0;i<NDIM;i++) // store indices of this cell 3063 2792 N[i] = LC->n[i]; 3064 2793 } else { 3065 DoeLog(1) && (eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl);2794 eLog() << Verbose(1) << "Point " << *a << " is not found in cell " << LC->index << "." << endl; 3066 2795 return; 3067 2796 } 3068 2797 // then go through the current and all neighbouring cells and check the contained points for possible candidates 3069 for (int i = 0; i < NDIM; i++) { 3070 Nlower[i] = ((N[i] - 1) >= 0) ? N[i] - 1 : 0; 3071 Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] - 1; 3072 } 3073 DoLog(0) && (Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl); 2798 for (int i=0;i<NDIM;i++) { 2799 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2800 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 2801 } 2802 Log() << Verbose(0) << "LC Intervals from [" << N[0] << "<->" << LC->N[0] << ", " << N[1] << "<->" << LC->N[1] << ", " << N[2] << "<->" << LC->N[2] << "] :" 2803 << " [" << Nlower[0] << "," << Nupper[0] << "], " << " [" << Nlower[1] << "," << Nupper[1] << "], " << " [" << Nlower[2] << "," << Nupper[2] << "], " << endl; 3074 2804 3075 2805 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 3076 2806 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 3077 2807 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 3078 const Linked Cell::LinkedNodes *List = LC->GetCurrentCell();2808 const LinkedNodes *List = LC->GetCurrentCell(); 3079 2809 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 3080 2810 if (List != NULL) { 3081 for (Linked Cell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {2811 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { 3082 2812 Candidate = (*Runner); 3083 2813 // check if we only have one unique point yet ... … … 3105 2835 norm = aCandidate.Norm(); 3106 2836 // second point shall have smallest angle with respect to Oben vector 3107 if (norm < RADIUS *2.) {2837 if (norm < RADIUS*2.) { 3108 2838 angle = AngleCheck.Angle(&Oben); 3109 2839 if (angle < Storage[0]) { 3110 2840 //Log() << Verbose(1) << "Old values of Storage: %lf %lf \n", Storage[0], Storage[1]); 3111 DoLog(1) && (Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n");2841 Log() << Verbose(1) << "Current candidate is " << *Candidate << ": Is a better candidate with distance " << norm << " and angle " << angle << " to oben " << Oben << ".\n"; 3112 2842 OptCandidate = Candidate; 3113 2843 Storage[0] = angle; … … 3124 2854 } 3125 2855 } else { 3126 DoLog(0) && (Log() << Verbose(0) << "Linked cell list is empty." << endl);2856 Log() << Verbose(0) << "Linked cell list is empty." << endl; 3127 2857 } 3128 2858 } 3129 } 3130 ; 2859 }; 2860 3131 2861 3132 2862 /** This recursive function finds a third point, to form a triangle with two given ones. … … 3156 2886 * @param OldSphereCenter center of sphere for base triangle, relative to center of BaseLine, giving null angle for the parameter circle 3157 2887 * @param CandidateLine CandidateForTesselation with the current base line and list of candidates and ShortestAngle 3158 * @param Third Pointthird point to avoid in search2888 * @param ThirdNode third point to avoid in search 3159 2889 * @param RADIUS radius of sphere 3160 2890 * @param *LC LinkedCell structure with neighbouring points 3161 2891 */ 3162 void Tesselation::FindThirdPointForTesselation( const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet * const ThirdPoint, const double RADIUS, const LinkedCell *LC) const3163 { 3164 Info FunctionInfo(__func__);3165 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers2892 void Tesselation::FindThirdPointForTesselation(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class TesselPoint * const ThirdNode, const double RADIUS, const LinkedCell *LC) const 2893 { 2894 Info FunctionInfo(__func__); 2895 Vector CircleCenter; // center of the circle, i.e. of the band of sphere's centers 3166 2896 Vector CirclePlaneNormal; // normal vector defining the plane this circle lives in 3167 2897 Vector SphereCenter; 3168 Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility3169 Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility3170 Vector NewNormalVector; // normal vector of the Candidate's triangle2898 Vector NewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, first possibility 2899 Vector OtherNewSphereCenter; // center of the sphere defined by the two points of BaseLine and the one of Candidate, second possibility 2900 Vector NewNormalVector; // normal vector of the Candidate's triangle 3171 2901 Vector helper, OptCandidateCenter, OtherOptCandidateCenter; 3172 2902 Vector RelativeOldSphereCenter; … … 3179 2909 TesselPoint *Candidate = NULL; 3180 2910 3181 DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl); 3182 3183 // copy old center 3184 CandidateLine.OldCenter.CopyVector(&OldSphereCenter); 3185 CandidateLine.ThirdPoint = ThirdPoint; 3186 CandidateLine.pointlist.clear(); 2911 Log() << Verbose(1) << "INFO: NormalVector of BaseTriangle is " << NormalVector << "." << endl; 3187 2912 3188 2913 // construct center of circle … … 3198 2923 RelativeOldSphereCenter.SubtractVector(&CircleCenter); 3199 2924 3200 // calculate squared radius TesselPoint *Third Point,f circle3201 radius = CirclePlaneNormal.NormSquared() /4.;3202 if (radius < RADIUS *RADIUS) {3203 CircleRadius = RADIUS *RADIUS - radius;2925 // calculate squared radius TesselPoint *ThirdNode,f circle 2926 radius = CirclePlaneNormal.NormSquared()/4.; 2927 if (radius < RADIUS*RADIUS) { 2928 CircleRadius = RADIUS*RADIUS - radius; 3204 2929 CirclePlaneNormal.Normalize(); 3205 DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);2930 Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 3206 2931 3207 2932 // test whether old center is on the band's plane 3208 2933 if (fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 3209 DoeLog(1) && (eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl);2934 eLog() << Verbose(1) << "Something's very wrong here: RelativeOldSphereCenter is not on the band's plane as desired by " << fabs(RelativeOldSphereCenter.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 3210 2935 RelativeOldSphereCenter.ProjectOntoPlane(&CirclePlaneNormal); 3211 2936 } 3212 2937 radius = RelativeOldSphereCenter.NormSquared(); 3213 2938 if (fabs(radius - CircleRadius) < HULLEPSILON) { 3214 DoLog(1) && (Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl);2939 Log() << Verbose(1) << "INFO: RelativeOldSphereCenter is at " << RelativeOldSphereCenter << "." << endl; 3215 2940 3216 2941 // check SearchDirection 3217 DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);3218 if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way!3219 DoeLog(1) && (eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl);2942 Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2943 if (fabs(RelativeOldSphereCenter.ScalarProduct(&SearchDirection)) > HULLEPSILON) { // rotated the wrong way! 2944 eLog() << Verbose(1) << "SearchDirection and RelativeOldSphereCenter are not orthogonal!" << endl; 3220 2945 } 3221 2946 3222 2947 // get cell for the starting point 3223 2948 if (LC->SetIndexToVector(&CircleCenter)) { 3224 for (int i = 0; i < NDIM;i++) // store indices of this cell3225 N[i] = LC->n[i];2949 for(int i=0;i<NDIM;i++) // store indices of this cell 2950 N[i] = LC->n[i]; 3226 2951 //Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 3227 2952 } else { 3228 DoeLog(1) && (eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl);2953 eLog() << Verbose(1) << "Vector " << CircleCenter << " is outside of LinkedCell's bounding box." << endl; 3229 2954 return; 3230 2955 } 3231 2956 // then go through the current and all neighbouring cells and check the contained points for possible candidates 3232 2957 //Log() << Verbose(1) << "LC Intervals:"; 3233 for (int i = 0; i < NDIM;i++) {3234 Nlower[i] = ((N[i] - 1) >= 0) ? N[i] -1 : 0;3235 Nupper[i] = ((N[i] + 1) < LC->N[i]) ? N[i] + 1 : LC->N[i] -1;2958 for (int i=0;i<NDIM;i++) { 2959 Nlower[i] = ((N[i]-1) >= 0) ? N[i]-1 : 0; 2960 Nupper[i] = ((N[i]+1) < LC->N[i]) ? N[i]+1 : LC->N[i]-1; 3236 2961 //Log() << Verbose(0) << " [" << Nlower[i] << "," << Nupper[i] << "] "; 3237 2962 } … … 3240 2965 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 3241 2966 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 3242 const Linked Cell::LinkedNodes *List = LC->GetCurrentCell();2967 const LinkedNodes *List = LC->GetCurrentCell(); 3243 2968 //Log() << Verbose(1) << "Current cell is " << LC->n[0] << ", " << LC->n[1] << ", " << LC->n[2] << " with No. " << LC->index << "." << endl; 3244 2969 if (List != NULL) { 3245 for (Linked Cell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {2970 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { 3246 2971 Candidate = (*Runner); 3247 2972 3248 2973 // check for three unique points 3249 DoLog(2) && (Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl);3250 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){2974 Log() << Verbose(2) << "INFO: Current Candidate is " << *Candidate << " for BaseLine " << *CandidateLine.BaseLine << " with OldSphereCenter " << OldSphereCenter << "." << endl; 2975 if ((Candidate != CandidateLine.BaseLine->endpoints[0]->node) && (Candidate != CandidateLine.BaseLine->endpoints[1]->node) ){ 3251 2976 3252 2977 // find center on the plane 3253 2978 GetCenterofCircumcircle(&NewPlaneCenter, *CandidateLine.BaseLine->endpoints[0]->node->node, *CandidateLine.BaseLine->endpoints[1]->node->node, *Candidate->node); 3254 DoLog(1) && (Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl); 3255 3256 if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node) && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON)) { 3257 DoLog(1) && (Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl); 2979 Log() << Verbose(1) << "INFO: NewPlaneCenter is " << NewPlaneCenter << "." << endl; 2980 2981 if (NewNormalVector.MakeNormalVector(CandidateLine.BaseLine->endpoints[0]->node->node, CandidateLine.BaseLine->endpoints[1]->node->node, Candidate->node) 2982 && (fabs(NewNormalVector.NormSquared()) > HULLEPSILON) 2983 ) { 2984 Log() << Verbose(1) << "INFO: NewNormalVector is " << NewNormalVector << "." << endl; 3258 2985 radius = CandidateLine.BaseLine->endpoints[0]->node->node->DistanceSquared(&NewPlaneCenter); 3259 DoLog(1) && (Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl);3260 DoLog(1) && (Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl);3261 DoLog(1) && (Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl);3262 if (radius < RADIUS *RADIUS) {2986 Log() << Verbose(1) << "INFO: CircleCenter is at " << CircleCenter << ", CirclePlaneNormal is " << CirclePlaneNormal << " with circle radius " << sqrt(CircleRadius) << "." << endl; 2987 Log() << Verbose(1) << "INFO: SearchDirection is " << SearchDirection << "." << endl; 2988 Log() << Verbose(1) << "INFO: Radius of CircumCenterCircle is " << radius << "." << endl; 2989 if (radius < RADIUS*RADIUS) { 3263 2990 otherradius = CandidateLine.BaseLine->endpoints[1]->node->node->DistanceSquared(&NewPlaneCenter); 3264 if (fabs(radius - otherradius) < HULLEPSILON) { 3265 // construct both new centers 3266 NewSphereCenter.CopyVector(&NewPlaneCenter); 3267 OtherNewSphereCenter.CopyVector(&NewPlaneCenter); 3268 helper.CopyVector(&NewNormalVector); 3269 helper.Scale(sqrt(RADIUS * RADIUS - radius)); 3270 DoLog(2) && (Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl); 3271 NewSphereCenter.AddVector(&helper); 3272 DoLog(2) && (Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl); 3273 // OtherNewSphereCenter is created by the same vector just in the other direction 3274 helper.Scale(-1.); 3275 OtherNewSphereCenter.AddVector(&helper); 3276 DoLog(2) && (Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl); 3277 3278 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 3279 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 3280 if ((ThirdPoint != NULL) && (Candidate == ThirdPoint->node)) { // in that case only the other circlecenter is valid 3281 if (OldSphereCenter.DistanceSquared(&NewSphereCenter) < OldSphereCenter.DistanceSquared(&OtherNewSphereCenter)) 3282 alpha = Otheralpha; 3283 } else 3284 alpha = min(alpha, Otheralpha); 3285 3286 // if there is a better candidate, drop the current list and add the new candidate 3287 // otherwise ignore the new candidate and keep the list 3288 if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) { 3289 if (fabs(alpha - Otheralpha) > MYEPSILON) { 3290 CandidateLine.OptCenter.CopyVector(&NewSphereCenter); 3291 CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter); 3292 } else { 3293 CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter); 3294 CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter); 3295 } 3296 // if there is an equal candidate, add it to the list without clearing the list 3297 if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) { 3298 CandidateLine.pointlist.push_back(Candidate); 3299 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl); 3300 } else { 3301 // remove all candidates from the list and then the list itself 3302 CandidateLine.pointlist.clear(); 3303 CandidateLine.pointlist.push_back(Candidate); 3304 DoLog(0) && (Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl); 3305 } 3306 CandidateLine.ShortestAngle = alpha; 3307 DoLog(0) && (Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl); 2991 if (fabs(radius - otherradius) > HULLEPSILON) { 2992 eLog() << Verbose(1) << "Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius-otherradius) << endl; 2993 } 2994 // construct both new centers 2995 NewSphereCenter.CopyVector(&NewPlaneCenter); 2996 OtherNewSphereCenter.CopyVector(&NewPlaneCenter); 2997 helper.CopyVector(&NewNormalVector); 2998 helper.Scale(sqrt(RADIUS*RADIUS - radius)); 2999 Log() << Verbose(2) << "INFO: Distance of NewPlaneCenter " << NewPlaneCenter << " to either NewSphereCenter is " << helper.Norm() << " of vector " << helper << " with sphere radius " << RADIUS << "." << endl; 3000 NewSphereCenter.AddVector(&helper); 3001 Log() << Verbose(2) << "INFO: NewSphereCenter is at " << NewSphereCenter << "." << endl; 3002 // OtherNewSphereCenter is created by the same vector just in the other direction 3003 helper.Scale(-1.); 3004 OtherNewSphereCenter.AddVector(&helper); 3005 Log() << Verbose(2) << "INFO: OtherNewSphereCenter is at " << OtherNewSphereCenter << "." << endl; 3006 3007 alpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, NewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 3008 Otheralpha = GetPathLengthonCircumCircle(CircleCenter, CirclePlaneNormal, CircleRadius, OtherNewSphereCenter, OldSphereCenter, NormalVector, SearchDirection); 3009 alpha = min(alpha, Otheralpha); 3010 3011 // if there is a better candidate, drop the current list and add the new candidate 3012 // otherwise ignore the new candidate and keep the list 3013 if (CandidateLine.ShortestAngle > (alpha - HULLEPSILON)) { 3014 if (fabs(alpha - Otheralpha) > MYEPSILON) { 3015 CandidateLine.OptCenter.CopyVector(&NewSphereCenter); 3016 CandidateLine.OtherOptCenter.CopyVector(&OtherNewSphereCenter); 3308 3017 } else { 3309 if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) { 3310 DoLog(1) && (Log() << Verbose(1) << "REJECT: Old candidate " << *(*CandidateLine.pointlist.begin()) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl); 3311 } else { 3312 DoLog(1) && (Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl); 3313 } 3018 CandidateLine.OptCenter.CopyVector(&OtherNewSphereCenter); 3019 CandidateLine.OtherOptCenter.CopyVector(&NewSphereCenter); 3314 3020 } 3021 // if there is an equal candidate, add it to the list without clearing the list 3022 if ((CandidateLine.ShortestAngle - HULLEPSILON) < alpha) { 3023 CandidateLine.pointlist.push_back(Candidate); 3024 Log() << Verbose(0) << "ACCEPT: We have found an equally good candidate: " << *(Candidate) << " with " 3025 << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl; 3026 } else { 3027 // remove all candidates from the list and then the list itself 3028 CandidateLine.pointlist.clear(); 3029 CandidateLine.pointlist.push_back(Candidate); 3030 Log() << Verbose(0) << "ACCEPT: We have found a better candidate: " << *(Candidate) << " with " 3031 << alpha << " and circumsphere's center at " << CandidateLine.OptCenter << "." << endl; 3032 } 3033 CandidateLine.ShortestAngle = alpha; 3034 Log() << Verbose(0) << "INFO: There are " << CandidateLine.pointlist.size() << " candidates in the list now." << endl; 3315 3035 } else { 3316 DoLog(1) && (Log() << Verbose(1) << "REJECT: Distance to center of circumcircle is not the same from each corner of the triangle: " << fabs(radius - otherradius) << endl); 3036 if ((Candidate != NULL) && (CandidateLine.pointlist.begin() != CandidateLine.pointlist.end())) { 3037 Log() << Verbose(1) << "REJECT: Old candidate " << *(Candidate) << " with " << CandidateLine.ShortestAngle << " is better than new one " << *Candidate << " with " << alpha << " ." << endl; 3038 } else { 3039 Log() << Verbose(1) << "REJECT: Candidate " << *Candidate << " with " << alpha << " was rejected." << endl; 3040 } 3317 3041 } 3318 3042 } else { 3319 DoLog(1) && (Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl);3043 Log() << Verbose(1) << "REJECT: NewSphereCenter " << NewSphereCenter << " for " << *Candidate << " is too far away: " << radius << "." << endl; 3320 3044 } 3321 3045 } else { 3322 DoLog(1) && (Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl);3046 Log() << Verbose(1) << "REJECT: Three points from " << *CandidateLine.BaseLine << " and Candidate " << *Candidate << " are linear-dependent." << endl; 3323 3047 } 3324 3048 } else { 3325 if (Third Point!= NULL) {3326 DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdPoint << " contains Candidate " << *Candidate << "." << endl);3049 if (ThirdNode != NULL) { 3050 Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " and " << *ThirdNode << " contains Candidate " << *Candidate << "." << endl; 3327 3051 } else { 3328 DoLog(1) && (Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl);3052 Log() << Verbose(1) << "REJECT: Base triangle " << *CandidateLine.BaseLine << " contains Candidate " << *Candidate << "." << endl; 3329 3053 } 3330 3054 } … … 3333 3057 } 3334 3058 } else { 3335 DoeLog(1) && (eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl);3059 eLog() << Verbose(1) << "The projected center of the old sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 3336 3060 } 3337 3061 } else { 3338 if (Third Point!= NULL)3339 DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdPoint << " is too big!" << endl);3062 if (ThirdNode != NULL) 3063 Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " and third node " << *ThirdNode << " is too big!" << endl; 3340 3064 else 3341 DoLog(1) && (Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl);3342 } 3343 3344 DoLog(1) && (Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl);3065 Log() << Verbose(1) << "Circumcircle for base line " << *CandidateLine.BaseLine << " is too big!" << endl; 3066 } 3067 3068 Log() << Verbose(1) << "INFO: Sorting candidate list ..." << endl; 3345 3069 if (CandidateLine.pointlist.size() > 1) { 3346 3070 CandidateLine.pointlist.unique(); 3347 3071 CandidateLine.pointlist.sort(); //SortCandidates); 3348 3072 } 3349 3350 if ((!CandidateLine.pointlist.empty()) && (!CandidateLine.CheckValidity(RADIUS, LC))) { 3351 DoeLog(0) && (eLog() << Verbose(0) << "There were other points contained in the rolling sphere as well!" << endl); 3352 performCriticalExit(); 3353 } 3354 } 3355 ; 3073 }; 3356 3074 3357 3075 /** Finds the endpoint two lines are sharing. … … 3362 3080 class BoundaryPointSet *Tesselation::GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const 3363 3081 { 3364 Info FunctionInfo(__func__);3082 Info FunctionInfo(__func__); 3365 3083 const BoundaryLineSet * lines[2] = { line1, line2 }; 3366 3084 class BoundaryPointSet *node = NULL; … … 3369 3087 for (int i = 0; i < 2; i++) 3370 3088 // for both lines 3371 for (int j = 0; j < 2; j++) { // for both endpoints 3372 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> (lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j])); 3373 if (!OrderTest.second) { // if insertion fails, we have common endpoint 3374 node = OrderTest.first->second; 3375 DoLog(1) && (Log() << Verbose(1) << "Common endpoint of lines " << *line1 << " and " << *line2 << " is: " << *node << "." << endl); 3376 j = 2; 3377 i = 2; 3378 break; 3089 for (int j = 0; j < 2; j++) 3090 { // for both endpoints 3091 OrderTest = OrderMap.insert(pair<int, class BoundaryPointSet *> ( 3092 lines[i]->endpoints[j]->Nr, lines[i]->endpoints[j])); 3093 if (!OrderTest.second) 3094 { // if insertion fails, we have common endpoint 3095 node = OrderTest.first->second; 3096 Log() << Verbose(1) << "Common endpoint of lines " << *line1 3097 << " and " << *line2 << " is: " << *node << "." << endl; 3098 j = 2; 3099 i = 2; 3100 break; 3101 } 3379 3102 } 3380 }3381 3103 return node; 3382 } 3383 ; 3104 }; 3384 3105 3385 3106 /** Finds the boundary points that are closest to a given Vector \a *x. … … 3395 3116 3396 3117 if (LinesOnBoundary.empty()) { 3397 DoeLog(1) && (eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl);3118 eLog() << Verbose(1) << "There is no tesselation structure to compare the point with, please create one first." << endl; 3398 3119 return NULL; 3399 3120 } … … 3401 3122 // gather all points close to the desired one 3402 3123 LC->SetIndexToVector(x); // ignore status as we calculate bounds below sensibly 3403 for (int i = 0; i < NDIM;i++) // store indices of this cell3124 for(int i=0;i<NDIM;i++) // store indices of this cell 3404 3125 N[i] = LC->n[i]; 3405 DoLog(1) && (Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl); 3126 Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 3127 3406 3128 DistanceToPointMap * points = new DistanceToPointMap; 3407 3129 LC->GetNeighbourBounds(Nlower, Nupper); … … 3410 3132 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 3411 3133 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 3412 const Linked Cell::LinkedNodes *List = LC->GetCurrentCell();3134 const LinkedNodes *List = LC->GetCurrentCell(); 3413 3135 //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl; 3414 3136 if (List != NULL) { 3415 for (Linked Cell::LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {3137 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { 3416 3138 FindPoint = PointsOnBoundary.find((*Runner)->nr); 3417 3139 if (FindPoint != PointsOnBoundary.end()) { 3418 points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second));3419 DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl);3140 points->insert(DistanceToPointPair (FindPoint->second->node->node->DistanceSquared(x), FindPoint->second) ); 3141 Log() << Verbose(1) << "INFO: Putting " << *FindPoint->second << " into the list." << endl; 3420 3142 } 3421 3143 } 3422 3144 } else { 3423 DoeLog(1) && (eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl);3145 eLog() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << " is invalid!" << endl; 3424 3146 } 3425 3147 } … … 3427 3149 // check whether we found some points 3428 3150 if (points->empty()) { 3429 DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);3430 delete (points);3151 eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl; 3152 delete(points); 3431 3153 return NULL; 3432 3154 } 3433 3155 return points; 3434 } 3435 ; 3156 }; 3436 3157 3437 3158 /** Finds the boundary line that is closest to a given Vector \a *x. … … 3443 3164 { 3444 3165 Info FunctionInfo(__func__); 3166 3445 3167 // get closest points 3446 DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x, LC);3168 DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC); 3447 3169 if (points == NULL) { 3448 DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);3170 eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl; 3449 3171 return NULL; 3450 3172 } 3451 3173 3452 3174 // for each point, check its lines, remember closest 3453 DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl);3175 Log() << Verbose(1) << "Finding closest BoundaryLine to " << *x << " ... " << endl; 3454 3176 BoundaryLineSet *ClosestLine = NULL; 3455 3177 double MinDistance = -1.; … … 3479 3201 helper.SubtractVector(&Center); 3480 3202 const double lengthB = helper.ScalarProduct(&BaseLine); 3481 if (lengthB * lengthA < 0) {// if have different sign3203 if (lengthB*lengthA < 0) { // if have different sign 3482 3204 ClosestLine = LineRunner->second; 3483 3205 MinDistance = distance; 3484 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl);3206 Log() << Verbose(1) << "ACCEPT: New closest line is " << *ClosestLine << " with projected distance " << MinDistance << "." << endl; 3485 3207 } else { 3486 DoLog(1) && (Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl);3208 Log() << Verbose(1) << "REJECT: Intersection is outside of the line section: " << lengthA << " and " << lengthB << "." << endl; 3487 3209 } 3488 3210 } else { 3489 DoLog(1) && (Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl);3211 Log() << Verbose(1) << "REJECT: Point is too further away than present line: " << distance << " >> " << MinDistance << "." << endl; 3490 3212 } 3491 3213 } 3492 3214 } 3493 delete (points);3215 delete(points); 3494 3216 // check whether closest line is "too close" :), then it's inside 3495 3217 if (ClosestLine == NULL) { 3496 DoLog(0) && (Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl);3218 Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl; 3497 3219 return NULL; 3498 3220 } 3499 3221 return ClosestLine; 3500 } 3501 ; 3222 }; 3223 3502 3224 3503 3225 /** Finds the triangle that is closest to a given Vector \a *x. … … 3508 3230 TriangleList * Tesselation::FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const 3509 3231 { 3510 Info FunctionInfo(__func__); 3511 // get closest points 3512 DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x, LC); 3232 Info FunctionInfo(__func__); 3233 3234 // get closest points 3235 DistanceToPointMap * points = FindClosestBoundaryPointsToVector(x,LC); 3513 3236 if (points == NULL) { 3514 DoeLog(1) && (eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl);3237 eLog() << Verbose(1) << "There is no nearest point: too far away from the surface." << endl; 3515 3238 return NULL; 3516 3239 } 3517 3240 3518 3241 // for each point, check its lines, remember closest 3519 DoLog(1) && (Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl);3242 Log() << Verbose(1) << "Finding closest BoundaryTriangle to " << *x << " ... " << endl; 3520 3243 LineSet ClosestLines; 3521 3244 double MinDistance = 1e+16; … … 3539 3262 const double lengthEndB = BaseLineIntersection.NormSquared(); 3540 3263 3541 if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) { // intersection would be outside, take closer endpoint3264 if ((lengthEndA > lengthBase) || (lengthEndB > lengthBase) || ((lengthEndA < MYEPSILON) || (lengthEndB < MYEPSILON))) { // intersection would be outside, take closer endpoint 3542 3265 const double lengthEnd = Min(lengthEndA, lengthEndB); 3543 3266 if (lengthEnd - MinDistance < -MYEPSILON) { // new best line … … 3545 3268 ClosestLines.insert(LineRunner->second); 3546 3269 MinDistance = lengthEnd; 3547 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl);3548 } else if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate3270 Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[0]->node << " is closer with " << lengthEnd << "." << endl; 3271 } else if (fabs(lengthEnd - MinDistance) < MYEPSILON) { // additional best candidate 3549 3272 ClosestLines.insert(LineRunner->second); 3550 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl);3273 Log() << Verbose(1) << "ACCEPT: Line " << *LineRunner->second << " to endpoint " << *LineRunner->second->endpoints[1]->node << " is equally good with " << lengthEnd << "." << endl; 3551 3274 } else { // line is worse 3552 DoLog(1) && (Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl);3275 Log() << Verbose(1) << "REJECT: Line " << *LineRunner->second << " to either endpoints is further away than present closest line candidate: " << lengthEndA << ", " << lengthEndB << ", and distance is longer than baseline:" << lengthBase << "." << endl; 3553 3276 } 3554 3277 } else { // intersection is closer, calculate … … 3561 3284 const double distance = BaseLineIntersection.NormSquared(); 3562 3285 if (Center.NormSquared() > BaseLine.NormSquared()) { 3563 DoeLog(0) && (eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl);3286 eLog() << Verbose(0) << "Algorithmic error: In second case we have intersection outside of baseline!" << endl; 3564 3287 } 3565 3288 if ((ClosestLines.empty()) || (distance < MinDistance)) { 3566 3289 ClosestLines.insert(LineRunner->second); 3567 3290 MinDistance = distance; 3568 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl);3291 Log() << Verbose(1) << "ACCEPT: Intersection in between endpoints, new closest line " << *LineRunner->second << " is " << *ClosestLines.begin() << " with projected distance " << MinDistance << "." << endl; 3569 3292 } else { 3570 DoLog(2) && (Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl);3293 Log() << Verbose(2) << "REJECT: Point is further away from line " << *LineRunner->second << " than present closest line: " << distance << " >> " << MinDistance << "." << endl; 3571 3294 } 3572 3295 } 3573 3296 } 3574 3297 } 3575 delete (points);3298 delete(points); 3576 3299 3577 3300 // check whether closest line is "too close" :), then it's inside 3578 3301 if (ClosestLines.empty()) { 3579 DoLog(0) && (Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl);3302 Log() << Verbose(0) << "Is the only point, no one else is closeby." << endl; 3580 3303 return NULL; 3581 3304 } … … 3583 3306 for (LineSet::iterator LineRunner = ClosestLines.begin(); LineRunner != ClosestLines.end(); LineRunner++) 3584 3307 for (TriangleMap::iterator Runner = (*LineRunner)->triangles.begin(); Runner != (*LineRunner)->triangles.end(); Runner++) { 3585 candidates->push_back(Runner->second);3586 }3308 candidates->push_back(Runner->second); 3309 } 3587 3310 return candidates; 3588 } 3589 ; 3311 }; 3590 3312 3591 3313 /** Finds closest triangle to a point. … … 3593 3315 * \param *out output stream for debugging 3594 3316 * \param *x Vector to look from 3595 * \param &distance contains found distance on return3596 3317 * \return list of BoundaryTriangleSet of nearest triangles or NULL. 3597 3318 */ 3598 3319 class BoundaryTriangleSet * Tesselation::FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const 3599 3320 { 3600 Info FunctionInfo(__func__);3321 Info FunctionInfo(__func__); 3601 3322 class BoundaryTriangleSet *result = NULL; 3602 3323 TriangleList *triangles = FindClosestTrianglesToVector(x, LC); … … 3609 3330 3610 3331 // go through all and pick the one with the best alignment to x 3611 double MinAlignment = 2. *M_PI;3332 double MinAlignment = 2.*M_PI; 3612 3333 for (TriangleList::iterator Runner = triangles->begin(); Runner != triangles->end(); Runner++) { 3613 3334 (*Runner)->GetCenter(&Center); … … 3618 3339 result = *Runner; 3619 3340 MinAlignment = Alignment; 3620 DoLog(1) && (Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl);3341 Log() << Verbose(1) << "ACCEPT: Triangle " << *result << " is better aligned with " << MinAlignment << "." << endl; 3621 3342 } else { 3622 DoLog(1) && (Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl);3623 } 3624 } 3625 delete (triangles);3343 Log() << Verbose(1) << "REJECT: Triangle " << *result << " is worse aligned with " << MinAlignment << "." << endl; 3344 } 3345 } 3346 delete(triangles); 3626 3347 3627 3348 return result; 3628 } 3629 ; 3349 }; 3630 3350 3631 3351 /** Checks whether the provided Vector is within the Tesselation structure. … … 3638 3358 bool Tesselation::IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const 3639 3359 { 3640 Info FunctionInfo(__func__); 3641 TriangleIntersectionList Intersections(&Point, this, LC); 3642 3643 return Intersections.IsInside(); 3360 return (GetDistanceSquaredToSurface(Point, LC) < MYEPSILON); 3644 3361 } 3645 ;3646 3362 3647 3363 /** Returns the distance to the surface given by the tesselation. … … 3673 3389 3674 3390 if (triangle == NULL) {// is boundary point or only point in point cloud? 3675 DoLog(1) && (Log() << Verbose(1) << "No triangle given!" << endl);3391 Log() << Verbose(1) << "No triangle given!" << endl; 3676 3392 return -1.; 3677 3393 } else { 3678 DoLog(1) && (Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl);3394 Log() << Verbose(1) << "INFO: Closest triangle found is " << *triangle << " with normal vector " << triangle->NormalVector << "." << endl; 3679 3395 } 3680 3396 3681 3397 triangle->GetCenter(&Center); 3682 DoLog(2) && (Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl);3398 Log() << Verbose(2) << "INFO: Central point of the triangle is " << Center << "." << endl; 3683 3399 DistanceToCenter.CopyVector(&Center); 3684 3400 DistanceToCenter.SubtractVector(&Point); 3685 DoLog(2) && (Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl);3401 Log() << Verbose(2) << "INFO: Vector from point to test to center is " << DistanceToCenter << "." << endl; 3686 3402 3687 3403 // check whether we are on boundary … … 3692 3408 Center.SubtractVector(&triangle->NormalVector); // points towards MolCenter 3693 3409 DistanceToCenter.AddVector(&triangle->NormalVector); // points outside 3694 DoLog(1) && (Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl);3410 Log() << Verbose(1) << "INFO: Calling Intersection with " << Center << " and " << DistanceToCenter << "." << endl; 3695 3411 if (triangle->GetIntersectionInsideTriangle(&Center, &DistanceToCenter, &Intersection)) { 3696 DoLog(1) && (Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl);3412 Log() << Verbose(1) << Point << " is inner point: sufficiently close to boundary, " << Intersection << "." << endl; 3697 3413 return 0.; 3698 3414 } else { 3699 DoLog(1) && (Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl);3415 Log() << Verbose(1) << Point << " is NOT an inner point: on triangle plane but outside of triangle bounds." << endl; 3700 3416 return false; 3701 3417 } … … 3703 3419 // calculate smallest distance 3704 3420 distance = triangle->GetClosestPointInsideTriangle(&Point, &Intersection); 3705 DoLog(1) && (Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl);3421 Log() << Verbose(1) << "Closest point on triangle is " << Intersection << "." << endl; 3706 3422 3707 3423 // then check direction to boundary 3708 3424 if (DistanceToCenter.ScalarProduct(&triangle->NormalVector) > MYEPSILON) { 3709 DoLog(1) && (Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl);3425 Log() << Verbose(1) << Point << " is an inner point, " << distance << " below surface." << endl; 3710 3426 return -distance; 3711 3427 } else { 3712 DoLog(1) && (Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl);3428 Log() << Verbose(1) << Point << " is NOT an inner point, " << distance << " above surface." << endl; 3713 3429 return +distance; 3714 3430 } 3715 3431 } 3716 } 3717 ; 3718 3719 /** Calculates minimum distance from \a&Point to a tesselated surface. 3432 }; 3433 3434 /** Calculates distance to a tesselated surface. 3720 3435 * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle(). 3721 3436 * \param &Point point to calculate distance from … … 3723 3438 * \return distance squared to closest point on surface 3724 3439 */ 3725 double Tesselation::GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const 3726 { 3727 Info FunctionInfo(__func__); 3728 TriangleIntersectionList Intersections(&Point, this, LC); 3729 3730 return Intersections.GetSmallestDistance(); 3731 } 3732 ; 3733 3734 /** Calculates minimum distance from \a&Point to a tesselated surface. 3735 * Combines \sa FindClosestTrianglesToVector() and \sa GetDistanceSquaredToTriangle(). 3736 * \param &Point point to calculate distance from 3737 * \param *LC needed for finding closest points fast 3738 * \return distance squared to closest point on surface 3739 */ 3740 BoundaryTriangleSet * Tesselation::GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const 3741 { 3742 Info FunctionInfo(__func__); 3743 TriangleIntersectionList Intersections(&Point, this, LC); 3744 3745 return Intersections.GetClosestTriangle(); 3746 } 3747 ; 3440 double Tesselation::GetDistanceSquaredToSurface(const Vector &Point, const LinkedCell* const LC) const 3441 { 3442 BoundaryTriangleSet *triangle = FindClosestTriangleToVector(&Point, LC); 3443 const double distance = GetDistanceSquaredToTriangle(Point, triangle); 3444 return Min(distance, LC->RADIUS); 3445 }; 3748 3446 3749 3447 /** Gets all points connected to the provided point by triangulation lines. … … 3755 3453 TesselPointSet * Tesselation::GetAllConnectedPoints(const TesselPoint* const Point) const 3756 3454 { 3757 Info FunctionInfo(__func__);3758 TesselPointSet *connectedPoints = new TesselPointSet;3455 Info FunctionInfo(__func__); 3456 TesselPointSet *connectedPoints = new TesselPointSet; 3759 3457 class BoundaryPointSet *ReferencePoint = NULL; 3760 3458 TesselPoint* current; 3761 3459 bool takePoint = false; 3460 3762 3461 // find the respective boundary point 3763 3462 PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr); … … 3765 3464 ReferencePoint = PointRunner->second; 3766 3465 } else { 3767 DoeLog(2) && (eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);3466 eLog() << Verbose(2) << "GetAllConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl; 3768 3467 ReferencePoint = NULL; 3769 3468 } … … 3771 3470 // little trick so that we look just through lines connect to the BoundaryPoint 3772 3471 // OR fall-back to look through all lines if there is no such BoundaryPoint 3773 const LineMap *Lines; 3774 ; 3472 const LineMap *Lines;; 3775 3473 if (ReferencePoint != NULL) 3776 3474 Lines = &(ReferencePoint->lines); … … 3779 3477 LineMap::const_iterator findLines = Lines->begin(); 3780 3478 while (findLines != Lines->end()) { 3781 takePoint = false;3782 3783 if (findLines->second->endpoints[0]->Nr == Point->nr) {3784 takePoint = true;3785 current = findLines->second->endpoints[1]->node;3786 } else if (findLines->second->endpoints[1]->Nr == Point->nr) {3787 takePoint = true;3788 current = findLines->second->endpoints[0]->node;3789 }3790 3791 if (takePoint) {3792 DoLog(1) && (Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl);3793 connectedPoints->insert(current);3794 }3795 3796 findLines++;3479 takePoint = false; 3480 3481 if (findLines->second->endpoints[0]->Nr == Point->nr) { 3482 takePoint = true; 3483 current = findLines->second->endpoints[1]->node; 3484 } else if (findLines->second->endpoints[1]->Nr == Point->nr) { 3485 takePoint = true; 3486 current = findLines->second->endpoints[0]->node; 3487 } 3488 3489 if (takePoint) { 3490 Log() << Verbose(1) << "INFO: Endpoint " << *current << " of line " << *(findLines->second) << " is enlisted." << endl; 3491 connectedPoints->insert(current); 3492 } 3493 3494 findLines++; 3797 3495 } 3798 3496 3799 3497 if (connectedPoints->empty()) { // if have not found any points 3800 DoeLog(1) && (eLog() << Verbose(1) << "We have not found any connected points to " << *Point << "." << endl);3498 eLog() << Verbose(1) << "We have not found any connected points to " << *Point<< "." << endl; 3801 3499 return NULL; 3802 3500 } 3803 3501 3804 3502 return connectedPoints; 3805 } 3806 ; 3503 }; 3504 3807 3505 3808 3506 /** Gets all points connected to the provided point by triangulation lines, ordered such that we have the circle round the point. … … 3820 3518 TesselPointList * Tesselation::GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const 3821 3519 { 3822 Info FunctionInfo(__func__);3520 Info FunctionInfo(__func__); 3823 3521 map<double, TesselPoint*> anglesOfPoints; 3824 3522 TesselPointList *connectedCircle = new TesselPointList; … … 3827 3525 Vector OrthogonalVector; 3828 3526 Vector helper; 3829 const TesselPoint * const TrianglePoints[3] = { Point, NULL, NULL};3527 const TesselPoint * const TrianglePoints[3] = {Point, NULL, NULL}; 3830 3528 TriangleList *triangles = NULL; 3831 3529 3832 3530 if (SetOfNeighbours == NULL) { 3833 DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);3834 delete (connectedCircle);3531 eLog() << Verbose(2) << "Could not find any connected points!" << endl; 3532 delete(connectedCircle); 3835 3533 return NULL; 3836 3534 } … … 3842 3540 PlaneNormal.AddVector(&(*Runner)->NormalVector); 3843 3541 } else { 3844 DoeLog(0) && (eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl);3542 eLog() << Verbose(0) << "Could not find any triangles for point " << *Point << "." << endl; 3845 3543 performCriticalExit(); 3846 3544 } 3847 PlaneNormal.Scale(1.0 /triangles->size());3848 DoLog(1) && (Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl);3545 PlaneNormal.Scale(1.0/triangles->size()); 3546 Log() << Verbose(1) << "INFO: Calculated PlaneNormal of all circle points is " << PlaneNormal << "." << endl; 3849 3547 PlaneNormal.Normalize(); 3850 3548 … … 3855 3553 AngleZero.ProjectOntoPlane(&PlaneNormal); 3856 3554 } 3857 if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {3858 DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);3555 if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) { 3556 Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl; 3859 3557 AngleZero.CopyVector((*SetOfNeighbours->begin())->node); 3860 3558 AngleZero.SubtractVector(Point->node); 3861 3559 AngleZero.ProjectOntoPlane(&PlaneNormal); 3862 3560 if (AngleZero.NormSquared() < MYEPSILON) { 3863 DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);3561 eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl; 3864 3562 performCriticalExit(); 3865 3563 } 3866 3564 } 3867 DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);3565 Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl; 3868 3566 if (AngleZero.NormSquared() > MYEPSILON) 3869 3567 OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero); 3870 3568 else 3871 3569 OrthogonalVector.MakeNormalVector(&PlaneNormal); 3872 DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);3570 Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl; 3873 3571 3874 3572 // go through all connected points and calculate angle … … 3878 3576 helper.ProjectOntoPlane(&PlaneNormal); 3879 3577 double angle = GetAngle(helper, AngleZero, OrthogonalVector); 3880 DoLog(0) && (Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl);3881 anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));3882 } 3883 3884 for (map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {3578 Log() << Verbose(0) << "INFO: Calculated angle is " << angle << " for point " << **listRunner << "." << endl; 3579 anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner))); 3580 } 3581 3582 for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) { 3885 3583 connectedCircle->push_back(AngleRunner->second); 3886 3584 } … … 3912 3610 3913 3611 if (SetOfNeighbours == NULL) { 3914 DoeLog(2) && (eLog() << Verbose(2) << "Could not find any connected points!" << endl);3915 delete (connectedCircle);3612 eLog() << Verbose(2) << "Could not find any connected points!" << endl; 3613 delete(connectedCircle); 3916 3614 return NULL; 3917 3615 } … … 3924 3622 } 3925 3623 3926 DoLog(1) && (Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl);3624 Log() << Verbose(1) << "INFO: Point is " << *Point << " and Reference is " << *Reference << "." << endl; 3927 3625 // calculate central point 3626 3928 3627 TesselPointSet::const_iterator TesselA = SetOfNeighbours->begin(); 3929 3628 TesselPointSet::const_iterator TesselB = SetOfNeighbours->begin(); … … 3935 3634 while (TesselC != SetOfNeighbours->end()) { 3936 3635 helper.MakeNormalVector((*TesselA)->node, (*TesselB)->node, (*TesselC)->node); 3937 DoLog(0) && (Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl);3636 Log() << Verbose(0) << "Making normal vector out of " << *(*TesselA) << ", " << *(*TesselB) << " and " << *(*TesselC) << ":" << helper << endl; 3938 3637 counter++; 3939 3638 TesselA++; … … 3944 3643 //Log() << Verbose(0) << "Summed vectors " << center << "; number of points " << connectedPoints.size() 3945 3644 // << "; scale factor " << counter; 3946 PlaneNormal.Scale(1.0 / (double)counter);3947 // Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl;3948 //3949 // // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points3950 // PlaneNormal.CopyVector(Point->node);3951 // PlaneNormal.SubtractVector(¢er);3952 // PlaneNormal.Normalize();3953 DoLog(1) && (Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl);3645 PlaneNormal.Scale(1.0/(double)counter); 3646 // Log() << Verbose(1) << "INFO: Calculated center of all circle points is " << center << "." << endl; 3647 // 3648 // // projection plane of the circle is at the closes Point and normal is pointing away from center of all circle points 3649 // PlaneNormal.CopyVector(Point->node); 3650 // PlaneNormal.SubtractVector(¢er); 3651 // PlaneNormal.Normalize(); 3652 Log() << Verbose(1) << "INFO: Calculated plane normal of circle is " << PlaneNormal << "." << endl; 3954 3653 3955 3654 // construct one orthogonal vector … … 3959 3658 AngleZero.ProjectOntoPlane(&PlaneNormal); 3960 3659 } 3961 if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) {3962 DoLog(1) && (Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl);3660 if ((Reference == NULL) || (AngleZero.NormSquared() < MYEPSILON )) { 3661 Log() << Verbose(1) << "Using alternatively " << *(*SetOfNeighbours->begin())->node << " as angle 0 referencer." << endl; 3963 3662 AngleZero.CopyVector((*SetOfNeighbours->begin())->node); 3964 3663 AngleZero.SubtractVector(Point->node); 3965 3664 AngleZero.ProjectOntoPlane(&PlaneNormal); 3966 3665 if (AngleZero.NormSquared() < MYEPSILON) { 3967 DoeLog(0) && (eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl);3666 eLog() << Verbose(0) << "CRITIAL: AngleZero is 0 even with alternative reference. The algorithm has to be changed here!" << endl; 3968 3667 performCriticalExit(); 3969 3668 } 3970 3669 } 3971 DoLog(1) && (Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl);3670 Log() << Verbose(1) << "INFO: Reference vector on this plane representing angle 0 is " << AngleZero << "." << endl; 3972 3671 if (AngleZero.NormSquared() > MYEPSILON) 3973 3672 OrthogonalVector.MakeNormalVector(&PlaneNormal, &AngleZero); 3974 3673 else 3975 3674 OrthogonalVector.MakeNormalVector(&PlaneNormal); 3976 DoLog(1) && (Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl);3675 Log() << Verbose(1) << "INFO: OrthogonalVector on plane is " << OrthogonalVector << "." << endl; 3977 3676 3978 3677 // go through all connected points and calculate angle 3979 pair <map<double, TesselPoint*>::iterator, bool> InserterTest;3678 pair <map<double, TesselPoint*>::iterator, bool > InserterTest; 3980 3679 for (TesselPointSet::iterator listRunner = SetOfNeighbours->begin(); listRunner != SetOfNeighbours->end(); listRunner++) { 3981 3680 helper.CopyVector((*listRunner)->node); … … 3984 3683 double angle = GetAngle(helper, AngleZero, OrthogonalVector); 3985 3684 if (angle > M_PI) // the correction is of no use here (and not desired) 3986 angle = 2. *M_PI - angle;3987 DoLog(0) && (Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl);3988 InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*> (angle, (*listRunner)));3685 angle = 2.*M_PI - angle; 3686 Log() << Verbose(0) << "INFO: Calculated angle between " << helper << " and " << AngleZero << " is " << angle << " for point " << **listRunner << "." << endl; 3687 InserterTest = anglesOfPoints.insert(pair<double, TesselPoint*>(angle, (*listRunner))); 3989 3688 if (!InserterTest.second) { 3990 DoeLog(0) && (eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl);3689 eLog() << Verbose(0) << "GetCircleOfSetOfPoints() got two atoms with same angle: " << *((InserterTest.first)->second) << " and " << (*listRunner) << endl; 3991 3690 performCriticalExit(); 3992 3691 } 3993 3692 } 3994 3693 3995 for (map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) {3694 for(map<double, TesselPoint*>::iterator AngleRunner = anglesOfPoints.begin(); AngleRunner != anglesOfPoints.end(); AngleRunner++) { 3996 3695 connectedCircle->push_back(AngleRunner->second); 3997 3696 } … … 4008 3707 ListOfTesselPointList * Tesselation::GetPathsOfConnectedPoints(const TesselPoint* const Point) const 4009 3708 { 4010 Info FunctionInfo(__func__);3709 Info FunctionInfo(__func__); 4011 3710 map<double, TesselPoint*> anglesOfPoints; 4012 list< TesselPointList *> *ListOfPaths = new list<TesselPointList *>;3711 list< TesselPointList *> *ListOfPaths = new list< TesselPointList *>; 4013 3712 TesselPointList *connectedPath = NULL; 4014 3713 Vector center; … … 4022 3721 class BoundaryLineSet *CurrentLine = NULL; 4023 3722 class BoundaryLineSet *StartLine = NULL; 3723 4024 3724 // find the respective boundary point 4025 3725 PointMap::const_iterator PointRunner = PointsOnBoundary.find(Point->nr); … … 4027 3727 ReferencePoint = PointRunner->second; 4028 3728 } else { 4029 DoeLog(1) && (eLog() << Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl);3729 eLog() << Verbose(1) << "GetPathOfConnectedPoints() could not find the BoundaryPoint belonging to " << *Point << "." << endl; 4030 3730 return NULL; 4031 3731 } 4032 3732 4033 map <class BoundaryLineSet *, bool> TouchedLine;4034 map <class BoundaryTriangleSet *, bool> TouchedTriangle;4035 map <class BoundaryLineSet *, bool>::iterator LineRunner;4036 map <class BoundaryTriangleSet *, bool>::iterator TriangleRunner;3733 map <class BoundaryLineSet *, bool> TouchedLine; 3734 map <class BoundaryTriangleSet *, bool> TouchedTriangle; 3735 map <class BoundaryLineSet *, bool>::iterator LineRunner; 3736 map <class BoundaryTriangleSet *, bool>::iterator TriangleRunner; 4037 3737 for (LineMap::iterator Runner = ReferencePoint->lines.begin(); Runner != ReferencePoint->lines.end(); Runner++) { 4038 TouchedLine.insert( pair<class BoundaryLineSet *, bool> (Runner->second, false));3738 TouchedLine.insert( pair <class BoundaryLineSet *, bool>(Runner->second, false) ); 4039 3739 for (TriangleMap::iterator Sprinter = Runner->second->triangles.begin(); Sprinter != Runner->second->triangles.end(); Sprinter++) 4040 TouchedTriangle.insert( pair<class BoundaryTriangleSet *, bool> (Sprinter->second, false));3740 TouchedTriangle.insert( pair <class BoundaryTriangleSet *, bool>(Sprinter->second, false) ); 4041 3741 } 4042 3742 if (!ReferencePoint->lines.empty()) { … … 4044 3744 LineRunner = TouchedLine.find(runner->second); 4045 3745 if (LineRunner == TouchedLine.end()) { 4046 DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl);3746 eLog() << Verbose(1) << "I could not find " << *runner->second << " in the touched list." << endl; 4047 3747 } else if (!LineRunner->second) { 4048 3748 LineRunner->second = true; … … 4052 3752 StartLine = CurrentLine; 4053 3753 CurrentPoint = CurrentLine->GetOtherEndpoint(ReferencePoint); 4054 DoLog(1) && (Log() << Verbose(1) << "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl);3754 Log() << Verbose(1)<< "INFO: Beginning path retrieval at " << *CurrentPoint << " of line " << *CurrentLine << "." << endl; 4055 3755 do { 4056 3756 // push current one 4057 DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl);3757 Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl; 4058 3758 connectedPath->push_back(CurrentPoint->node); 4059 3759 4060 3760 // find next triangle 4061 3761 for (TriangleMap::iterator Runner = CurrentLine->triangles.begin(); Runner != CurrentLine->triangles.end(); Runner++) { 4062 DoLog(1) && (Log() << Verbose(1) << "INFO: Inspecting triangle " << *Runner->second << "." << endl);3762 Log() << Verbose(1) << "INFO: Inspecting triangle " << *Runner->second << "." << endl; 4063 3763 if ((Runner->second != triangle)) { // look for first triangle not equal to old one 4064 3764 triangle = Runner->second; … … 4067 3767 if (!TriangleRunner->second) { 4068 3768 TriangleRunner->second = true; 4069 DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl);3769 Log() << Verbose(1) << "INFO: Connecting triangle is " << *triangle << "." << endl; 4070 3770 break; 4071 3771 } else { 4072 DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl);3772 Log() << Verbose(1) << "INFO: Skipping " << *triangle << ", as we have already visited it." << endl; 4073 3773 triangle = NULL; 4074 3774 } 4075 3775 } else { 4076 DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl);3776 eLog() << Verbose(1) << "I could not find " << *triangle << " in the touched list." << endl; 4077 3777 triangle = NULL; 4078 3778 } … … 4082 3782 break; 4083 3783 // find next line 4084 for (int i = 0; i < 3;i++) {3784 for (int i=0;i<3;i++) { 4085 3785 if ((triangle->lines[i] != CurrentLine) && (triangle->lines[i]->ContainsBoundaryPoint(ReferencePoint))) { // not the current line and still containing Point 4086 3786 CurrentLine = triangle->lines[i]; 4087 DoLog(1) && (Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl);3787 Log() << Verbose(1) << "INFO: Connecting line is " << *CurrentLine << "." << endl; 4088 3788 break; 4089 3789 } … … 4091 3791 LineRunner = TouchedLine.find(CurrentLine); 4092 3792 if (LineRunner == TouchedLine.end()) 4093 DoeLog(1) && (eLog() << Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl);3793 eLog() << Verbose(1) << "I could not find " << *CurrentLine << " in the touched list." << endl; 4094 3794 else 4095 3795 LineRunner->second = true; … … 4099 3799 } while (CurrentLine != StartLine); 4100 3800 // last point is missing, as it's on start line 4101 DoLog(1) && (Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl);3801 Log() << Verbose(1) << "INFO: Putting " << *CurrentPoint << " at end of path." << endl; 4102 3802 if (StartLine->GetOtherEndpoint(ReferencePoint)->node != connectedPath->back()) 4103 3803 connectedPath->push_back(StartLine->GetOtherEndpoint(ReferencePoint)->node); … … 4105 3805 ListOfPaths->push_back(connectedPath); 4106 3806 } else { 4107 DoLog(1) && (Log() << Verbose(1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl);3807 Log() << Verbose(1) << "INFO: Skipping " << *runner->second << ", as we have already visited it." << endl; 4108 3808 } 4109 3809 } 4110 3810 } else { 4111 DoeLog(1) && (eLog() << Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl);3811 eLog() << Verbose(1) << "There are no lines attached to " << *ReferencePoint << "." << endl; 4112 3812 } 4113 3813 … … 4123 3823 ListOfTesselPointList * Tesselation::GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const 4124 3824 { 4125 Info FunctionInfo(__func__);3825 Info FunctionInfo(__func__); 4126 3826 list<TesselPointList *> *ListofPaths = GetPathsOfConnectedPoints(Point); 4127 list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *> ;3827 list<TesselPointList *> *ListofClosedPaths = new list<TesselPointList *>; 4128 3828 TesselPointList *connectedPath = NULL; 4129 3829 TesselPointList *newPath = NULL; 4130 3830 int count = 0; 3831 3832 4131 3833 TesselPointList::iterator CircleRunner; 4132 3834 TesselPointList::iterator CircleStart; 4133 3835 4134 for (list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) {3836 for(list<TesselPointList *>::iterator ListRunner = ListofPaths->begin(); ListRunner != ListofPaths->end(); ListRunner++) { 4135 3837 connectedPath = *ListRunner; 4136 3838 4137 DoLog(1) && (Log() << Verbose(1) << "INFO: Current path is " << connectedPath << "." << endl);3839 Log() << Verbose(1) << "INFO: Current path is " << connectedPath << "." << endl; 4138 3840 4139 3841 // go through list, look for reappearance of starting Point and count 4140 3842 CircleStart = connectedPath->begin(); 3843 4141 3844 // go through list, look for reappearance of starting Point and create list 4142 3845 TesselPointList::iterator Marker = CircleStart; … … 4144 3847 if ((*CircleRunner == *CircleStart) && (CircleRunner != CircleStart)) { // is not the very first point 4145 3848 // we have a closed circle from Marker to new Marker 4146 DoLog(1) && (Log() << Verbose(1) << count + 1 << ". closed path consists of: ");3849 Log() << Verbose(1) << count+1 << ". closed path consists of: "; 4147 3850 newPath = new TesselPointList; 4148 3851 TesselPointList::iterator CircleSprinter = Marker; 4149 3852 for (; CircleSprinter != CircleRunner; CircleSprinter++) { 4150 3853 newPath->push_back(*CircleSprinter); 4151 DoLog(0) && (Log() << Verbose(0) << (**CircleSprinter) << " <-> ");3854 Log() << Verbose(0) << (**CircleSprinter) << " <-> "; 4152 3855 } 4153 DoLog(0) && (Log() << Verbose(0) << ".." << endl);3856 Log() << Verbose(0) << ".." << endl; 4154 3857 count++; 4155 3858 Marker = CircleRunner; … … 4160 3863 } 4161 3864 } 4162 DoLog(1) && (Log() << Verbose(1) << "INFO: " << count << " closed additional path(s) have been created." << endl);3865 Log() << Verbose(1) << "INFO: " << count << " closed additional path(s) have been created." << endl; 4163 3866 4164 3867 // delete list of paths … … 4166 3869 connectedPath = *(ListofPaths->begin()); 4167 3870 ListofPaths->remove(connectedPath); 4168 delete (connectedPath);4169 } 4170 delete (ListofPaths);3871 delete(connectedPath); 3872 } 3873 delete(ListofPaths); 4171 3874 4172 3875 // exit 4173 3876 return ListofClosedPaths; 4174 } 4175 ; 3877 }; 3878 4176 3879 4177 3880 /** Gets all belonging triangles for a given BoundaryPointSet. … … 4182 3885 TriangleSet *Tesselation::GetAllTriangles(const BoundaryPointSet * const Point) const 4183 3886 { 4184 Info FunctionInfo(__func__);4185 TriangleSet *connectedTriangles = new TriangleSet;3887 Info FunctionInfo(__func__); 3888 TriangleSet *connectedTriangles = new TriangleSet; 4186 3889 4187 3890 if (Point == NULL) { 4188 DoeLog(1) && (eLog() << Verbose(1) << "Point given is NULL." << endl);3891 eLog() << Verbose(1) << "Point given is NULL." << endl; 4189 3892 } else { 4190 3893 // go through its lines and insert all triangles 4191 3894 for (LineMap::const_iterator LineRunner = Point->lines.begin(); LineRunner != Point->lines.end(); LineRunner++) 4192 3895 for (TriangleMap::iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) { 4193 connectedTriangles->insert(TriangleRunner->second);4194 }3896 connectedTriangles->insert(TriangleRunner->second); 3897 } 4195 3898 } 4196 3899 4197 3900 return connectedTriangles; 4198 } 4199 ; 3901 }; 3902 4200 3903 4201 3904 /** Removes a boundary point from the envelope while keeping it closed. … … 4210 3913 * \return volume added to the volume inside the tesselated surface by the removal 4211 3914 */ 4212 double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point) 4213 { 3915 double Tesselation::RemovePointFromTesselatedSurface(class BoundaryPointSet *point) { 4214 3916 class BoundaryLineSet *line = NULL; 4215 3917 class BoundaryTriangleSet *triangle = NULL; … … 4219 3921 4220 3922 if (point == NULL) { 4221 DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl);3923 eLog() << Verbose(1) << "Cannot remove the point " << point << ", it's NULL!" << endl; 4222 3924 return 0.; 4223 3925 } else 4224 DoLog(0) && (Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl);3926 Log() << Verbose(0) << "Removing point " << *point << " from tesselated boundary ..." << endl; 4225 3927 4226 3928 // copy old location for the volume … … 4229 3931 // get list of connected points 4230 3932 if (point->lines.empty()) { 4231 DoeLog(1) && (eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl);3933 eLog() << Verbose(1) << "Cannot remove the point " << *point << ", it's connected to no lines!" << endl; 4232 3934 return 0.; 4233 3935 } … … 4238 3940 // gather all triangles 4239 3941 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) 4240 count +=LineRunner->second->triangles.size();3942 count+=LineRunner->second->triangles.size(); 4241 3943 TriangleMap Candidates; 4242 3944 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) { … … 4244 3946 for (TriangleMap::iterator TriangleRunner = line->triangles.begin(); TriangleRunner != line->triangles.end(); TriangleRunner++) { 4245 3947 triangle = TriangleRunner->second; 4246 Candidates.insert( TrianglePair(triangle->Nr, triangle));3948 Candidates.insert( TrianglePair (triangle->Nr, triangle) ); 4247 3949 } 4248 3950 } 4249 3951 4250 3952 // remove all triangles 4251 count =0;3953 count=0; 4252 3954 NormalVector.Zero(); 4253 3955 for (TriangleMap::iterator Runner = Candidates.begin(); Runner != Candidates.end(); Runner++) { 4254 DoLog(1) && (Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl);3956 Log() << Verbose(1) << "INFO: Removing triangle " << *(Runner->second) << "." << endl; 4255 3957 NormalVector.SubtractVector(&Runner->second->NormalVector); // has to point inward 4256 3958 RemoveTesselationTriangle(Runner->second); 4257 3959 count++; 4258 3960 } 4259 DoLog(1) && (Log() << Verbose(1) << count << " triangles were removed." << endl);3961 Log() << Verbose(1) << count << " triangles were removed." << endl; 4260 3962 4261 3963 list<TesselPointList *>::iterator ListAdvance = ListOfClosedPaths->begin(); … … 4266 3968 double smallestangle; 4267 3969 Vector Point, Reference, OrthogonalVector; 4268 if (count > 2) { // less than three triangles, then nothing will be created3970 if (count > 2) { // less than three triangles, then nothing will be created 4269 3971 class TesselPoint *TriangleCandidates[3]; 4270 3972 count = 0; 4271 for ( ; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) {// go through all closed paths3973 for ( ; ListRunner != ListOfClosedPaths->end(); ListRunner = ListAdvance) { // go through all closed paths 4272 3974 if (ListAdvance != ListOfClosedPaths->end()) 4273 3975 ListAdvance++; 4274 3976 4275 3977 connectedPath = *ListRunner; 3978 4276 3979 // re-create all triangles by going through connected points list 4277 3980 LineList NewLines; 4278 for (; !connectedPath->empty();) {3981 for (;!connectedPath->empty();) { 4279 3982 // search middle node with widest angle to next neighbours 4280 3983 EndNode = connectedPath->end(); 4281 3984 smallestangle = 0.; 4282 3985 for (MiddleNode = connectedPath->begin(); MiddleNode != connectedPath->end(); MiddleNode++) { 4283 DoLog(1) && (Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);3986 Log() << Verbose(1) << "INFO: MiddleNode is " << **MiddleNode << "." << endl; 4284 3987 // construct vectors to next and previous neighbour 4285 3988 StartNode = MiddleNode; … … 4302 4005 angle = GetAngle(Point, Reference, OrthogonalVector); 4303 4006 //if (angle < M_PI) // no wrong-sided triangles, please? 4304 if (fabs(angle - M_PI) < fabs(smallestangle - M_PI)) {// get straightest angle (i.e. construct those triangles with smallest area first)4305 smallestangle = angle;4306 EndNode = MiddleNode;4307 }4007 if(fabs(angle - M_PI) < fabs(smallestangle - M_PI)) { // get straightest angle (i.e. construct those triangles with smallest area first) 4008 smallestangle = angle; 4009 EndNode = MiddleNode; 4010 } 4308 4011 } 4309 4012 MiddleNode = EndNode; 4310 4013 if (MiddleNode == connectedPath->end()) { 4311 DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl);4014 eLog() << Verbose(0) << "CRITICAL: Could not find a smallest angle!" << endl; 4312 4015 performCriticalExit(); 4313 4016 } … … 4319 4022 if (EndNode == connectedPath->end()) 4320 4023 EndNode = connectedPath->begin(); 4321 DoLog(2) && (Log() << Verbose(2) << "INFO: StartNode is " << **StartNode << "." << endl);4322 DoLog(2) && (Log() << Verbose(2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl);4323 DoLog(2) && (Log() << Verbose(2) << "INFO: EndNode is " << **EndNode << "." << endl);4324 DoLog(1) && (Log() << Verbose(1) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl);4024 Log() << Verbose(2) << "INFO: StartNode is " << **StartNode << "." << endl; 4025 Log() << Verbose(2) << "INFO: MiddleNode is " << **MiddleNode << "." << endl; 4026 Log() << Verbose(2) << "INFO: EndNode is " << **EndNode << "." << endl; 4027 Log() << Verbose(1) << "INFO: Attempting to create triangle " << (*StartNode)->Name << ", " << (*MiddleNode)->Name << " and " << (*EndNode)->Name << "." << endl; 4325 4028 TriangleCandidates[0] = *StartNode; 4326 4029 TriangleCandidates[1] = *MiddleNode; … … 4328 4031 triangle = GetPresentTriangle(TriangleCandidates); 4329 4032 if (triangle != NULL) { 4330 DoeLog(0) && (eLog() << Verbose(0) << "New triangle already present, skipping!" << endl);4033 eLog() << Verbose(0) << "New triangle already present, skipping!" << endl; 4331 4034 StartNode++; 4332 4035 MiddleNode++; … … 4340 4043 continue; 4341 4044 } 4342 DoLog(3) && (Log() << Verbose(3) << "Adding new triangle points." << endl);4045 Log() << Verbose(3) << "Adding new triangle points."<< endl; 4343 4046 AddTesselationPoint(*StartNode, 0); 4344 4047 AddTesselationPoint(*MiddleNode, 1); 4345 4048 AddTesselationPoint(*EndNode, 2); 4346 DoLog(3) && (Log() << Verbose(3) << "Adding new triangle lines." << endl);4347 AddTesselationLine( NULL, NULL,TPS[0], TPS[1], 0);4348 AddTesselationLine( NULL, NULL,TPS[0], TPS[2], 1);4049 Log() << Verbose(3) << "Adding new triangle lines."<< endl; 4050 AddTesselationLine(TPS[0], TPS[1], 0); 4051 AddTesselationLine(TPS[0], TPS[2], 1); 4349 4052 NewLines.push_back(BLS[1]); 4350 AddTesselationLine( NULL, NULL,TPS[1], TPS[2], 2);4053 AddTesselationLine(TPS[1], TPS[2], 2); 4351 4054 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 4352 4055 BTS->GetNormalVector(NormalVector); … … 4359 4062 // prepare nodes for next triangle 4360 4063 StartNode = EndNode; 4361 DoLog(2) && (Log() << Verbose(2) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl);4064 Log() << Verbose(2) << "Removing " << **MiddleNode << " from closed path, remaining points: " << connectedPath->size() << "." << endl; 4362 4065 connectedPath->remove(*MiddleNode); // remove the middle node (it is surrounded by triangles) 4363 4066 if (connectedPath->size() == 2) { // we are done … … 4366 4069 break; 4367 4070 } else if (connectedPath->size() < 2) { // something's gone wrong! 4368 DoeLog(0) && (eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl);4071 eLog() << Verbose(0) << "CRITICAL: There are only two endpoints left!" << endl; 4369 4072 performCriticalExit(); 4370 4073 } else { … … 4386 4089 do { 4387 4090 maxgain = 0; 4388 for (LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) {4091 for(LineList::iterator Runner = NewLines.begin(); Runner != NewLines.end(); Runner++) { 4389 4092 tmp = PickFarthestofTwoBaselines(*Runner); 4390 4093 if (maxgain < tmp) { … … 4395 4098 if (maxgain != 0) { 4396 4099 volume += maxgain; 4397 DoLog(1) && (Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl);4100 Log() << Verbose(1) << "Flipping baseline with highest volume" << **Candidate << "." << endl; 4398 4101 OtherBase = FlipBaseline(*Candidate); 4399 4102 NewLines.erase(Candidate); … … 4404 4107 4405 4108 ListOfClosedPaths->remove(connectedPath); 4406 delete (connectedPath);4407 } 4408 DoLog(0) && (Log() << Verbose(0) << count << " triangles were created." << endl);4109 delete(connectedPath); 4110 } 4111 Log() << Verbose(0) << count << " triangles were created." << endl; 4409 4112 } else { 4410 4113 while (!ListOfClosedPaths->empty()) { … … 4412 4115 connectedPath = *ListRunner; 4413 4116 ListOfClosedPaths->remove(connectedPath); 4414 delete (connectedPath);4415 } 4416 DoLog(0) && (Log() << Verbose(0) << "No need to create any triangles." << endl);4417 } 4418 delete (ListOfClosedPaths);4419 4420 DoLog(0) && (Log() << Verbose(0) << "Removed volume is " << volume << "." << endl);4117 delete(connectedPath); 4118 } 4119 Log() << Verbose(0) << "No need to create any triangles." << endl; 4120 } 4121 delete(ListOfClosedPaths); 4122 4123 Log() << Verbose(0) << "Removed volume is " << volume << "." << endl; 4421 4124 4422 4125 return volume; 4423 } 4424 ; 4126 }; 4127 4128 4425 4129 4426 4130 /** … … 4434 4138 TriangleList *Tesselation::FindTriangles(const TesselPoint* const Points[3]) const 4435 4139 { 4436 Info FunctionInfo(__func__);4437 TriangleList *result = new TriangleList;4140 Info FunctionInfo(__func__); 4141 TriangleList *result = new TriangleList; 4438 4142 LineMap::const_iterator FindLine; 4439 4143 TriangleMap::const_iterator FindTriangle; … … 4459 4163 for (int i = 0; i < 3; i++) { 4460 4164 if (TrianglePoints[i] != NULL) { 4461 for (int j = i +1; j < 3; j++) {4165 for (int j = i+1; j < 3; j++) { 4462 4166 if (TrianglePoints[j] != NULL) { 4463 4167 for (FindLine = TrianglePoints[i]->lines.find(TrianglePoints[j]->node->nr); // is a multimap! 4464 (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr); FindLine++) { 4465 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) { 4168 (FindLine != TrianglePoints[i]->lines.end()) && (FindLine->first == TrianglePoints[j]->node->nr); 4169 FindLine++) { 4170 for (FindTriangle = FindLine->second->triangles.begin(); 4171 FindTriangle != FindLine->second->triangles.end(); 4172 FindTriangle++) { 4466 4173 if (FindTriangle->second->IsPresentTupel(TrianglePoints)) { 4467 4174 result->push_back(FindTriangle->second); … … 4478 4185 case 1: // copy all triangles of the respective line 4479 4186 { 4480 int i =0;4187 int i=0; 4481 4188 for (; i < 3; i++) 4482 4189 if (TrianglePoints[i] == NULL) 4483 4190 break; 4484 for (FindLine = TrianglePoints[(i + 1) % 3]->lines.find(TrianglePoints[(i + 2) % 3]->node->nr); // is a multimap! 4485 (FindLine != TrianglePoints[(i + 1) % 3]->lines.end()) && (FindLine->first == TrianglePoints[(i + 2) % 3]->node->nr); FindLine++) { 4486 for (FindTriangle = FindLine->second->triangles.begin(); FindTriangle != FindLine->second->triangles.end(); FindTriangle++) { 4191 for (FindLine = TrianglePoints[(i+1)%3]->lines.find(TrianglePoints[(i+2)%3]->node->nr); // is a multimap! 4192 (FindLine != TrianglePoints[(i+1)%3]->lines.end()) && (FindLine->first == TrianglePoints[(i+2)%3]->node->nr); 4193 FindLine++) { 4194 for (FindTriangle = FindLine->second->triangles.begin(); 4195 FindTriangle != FindLine->second->triangles.end(); 4196 FindTriangle++) { 4487 4197 if (FindTriangle->second->IsPresentTupel(TrianglePoints)) { 4488 4198 result->push_back(FindTriangle->second); … … 4494 4204 case 2: // copy all triangles of the respective point 4495 4205 { 4496 int i =0;4206 int i=0; 4497 4207 for (; i < 3; i++) 4498 4208 if (TrianglePoints[i] != NULL) … … 4512 4222 } 4513 4223 default: 4514 DoeLog(0) && (eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl);4224 eLog() << Verbose(0) << "Number of wildcards is greater than 3, cannot happen!" << endl; 4515 4225 performCriticalExit(); 4516 4226 break; … … 4520 4230 } 4521 4231 4522 struct BoundaryLineSetCompare 4523 { 4524 bool operator()(const BoundaryLineSet * const a, const BoundaryLineSet * const b) 4525 { 4232 struct BoundaryLineSetCompare { 4233 bool operator() (const BoundaryLineSet * const a, const BoundaryLineSet * const b) { 4526 4234 int lowerNra = -1; 4527 4235 int lowerNrb = -1; … … 4541 4249 else if (a->endpoints[lowerNra] > b->endpoints[lowerNrb]) 4542 4250 return false; 4543 else { // both lower-numbered endpoints are the same ...4544 if (a->endpoints[(lowerNra + 1) % 2] < b->endpoints[(lowerNrb + 1) %2])4545 return true;4546 else if (a->endpoints[(lowerNra + 1) % 2] > b->endpoints[(lowerNrb + 1) %2])4547 return false;4251 else { // both lower-numbered endpoints are the same ... 4252 if (a->endpoints[(lowerNra+1)%2] < b->endpoints[(lowerNrb+1)%2]) 4253 return true; 4254 else if (a->endpoints[(lowerNra+1)%2] > b->endpoints[(lowerNrb+1)%2]) 4255 return false; 4548 4256 } 4549 4257 return false; 4550 } 4551 ; 4258 }; 4552 4259 }; 4553 4260 … … 4562 4269 IndexToIndex * Tesselation::FindAllDegeneratedLines() 4563 4270 { 4564 Info FunctionInfo(__func__);4565 UniqueLines AllLines;4271 Info FunctionInfo(__func__); 4272 UniqueLines AllLines; 4566 4273 IndexToIndex * DegeneratedLines = new IndexToIndex; 4567 4274 4568 4275 // sanity check 4569 4276 if (LinesOnBoundary.empty()) { 4570 DoeLog(2) && (eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure.");4277 eLog() << Verbose(2) << "FindAllDegeneratedTriangles() was called without any tesselation structure."; 4571 4278 return DegeneratedLines; 4572 4279 } 4280 4573 4281 LineMap::iterator LineRunner1; 4574 pair< UniqueLines::iterator, bool> tester;4282 pair< UniqueLines::iterator, bool> tester; 4575 4283 for (LineRunner1 = LinesOnBoundary.begin(); LineRunner1 != LinesOnBoundary.end(); ++LineRunner1) { 4576 tester = AllLines.insert( LineRunner1->second);4284 tester = AllLines.insert( LineRunner1->second ); 4577 4285 if (!tester.second) { // found degenerated line 4578 DegeneratedLines->insert (pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr));4579 DegeneratedLines->insert (pair<int, int> ((*tester.first)->Nr, LineRunner1->second->Nr));4286 DegeneratedLines->insert ( pair<int, int> (LineRunner1->second->Nr, (*tester.first)->Nr) ); 4287 DegeneratedLines->insert ( pair<int, int> ((*tester.first)->Nr, LineRunner1->second->Nr) ); 4580 4288 } 4581 4289 } … … 4583 4291 AllLines.clear(); 4584 4292 4585 DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl);4293 Log() << Verbose(0) << "FindAllDegeneratedLines() found " << DegeneratedLines->size() << " lines." << endl; 4586 4294 IndexToIndex::iterator it; 4587 4295 for (it = DegeneratedLines->begin(); it != DegeneratedLines->end(); it++) { … … 4589 4297 const LineMap::const_iterator Line2 = LinesOnBoundary.find((*it).second); 4590 4298 if (Line1 != LinesOnBoundary.end() && Line2 != LinesOnBoundary.end()) 4591 DoLog(0) && (Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl);4299 Log() << Verbose(0) << *Line1->second << " => " << *Line2->second << endl; 4592 4300 else 4593 DoeLog(1) && (eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl);4301 eLog() << Verbose(1) << "Either " << (*it).first << " or " << (*it).second << " are not in LinesOnBoundary!" << endl; 4594 4302 } 4595 4303 … … 4605 4313 IndexToIndex * Tesselation::FindAllDegeneratedTriangles() 4606 4314 { 4607 Info FunctionInfo(__func__);4315 Info FunctionInfo(__func__); 4608 4316 IndexToIndex * DegeneratedLines = FindAllDegeneratedLines(); 4609 4317 IndexToIndex * DegeneratedTriangles = new IndexToIndex; 4318 4610 4319 TriangleMap::iterator TriangleRunner1, TriangleRunner2; 4611 4320 LineMap::iterator Liner; … … 4622 4331 for (TriangleRunner1 = line1->triangles.begin(); TriangleRunner1 != line1->triangles.end(); ++TriangleRunner1) { 4623 4332 for (TriangleRunner2 = line2->triangles.begin(); TriangleRunner2 != line2->triangles.end(); ++TriangleRunner2) { 4624 if ((TriangleRunner1->second != TriangleRunner2->second) && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) { 4625 DegeneratedTriangles->insert(pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr)); 4626 DegeneratedTriangles->insert(pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr)); 4333 if ((TriangleRunner1->second != TriangleRunner2->second) 4334 && (TriangleRunner1->second->IsPresentTupel(TriangleRunner2->second))) { 4335 DegeneratedTriangles->insert( pair<int, int> (TriangleRunner1->second->Nr, TriangleRunner2->second->Nr) ); 4336 DegeneratedTriangles->insert( pair<int, int> (TriangleRunner2->second->Nr, TriangleRunner1->second->Nr) ); 4627 4337 } 4628 4338 } 4629 4339 } 4630 4340 } 4631 delete (DegeneratedLines);4632 4633 DoLog(0) && (Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl);4341 delete(DegeneratedLines); 4342 4343 Log() << Verbose(0) << "FindAllDegeneratedTriangles() found " << DegeneratedTriangles->size() << " triangles:" << endl; 4634 4344 IndexToIndex::iterator it; 4635 4345 for (it = DegeneratedTriangles->begin(); it != DegeneratedTriangles->end(); it++) 4636 DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl);4346 Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl; 4637 4347 4638 4348 return DegeneratedTriangles; … … 4645 4355 void Tesselation::RemoveDegeneratedTriangles() 4646 4356 { 4647 Info FunctionInfo(__func__);4357 Info FunctionInfo(__func__); 4648 4358 IndexToIndex * DegeneratedTriangles = FindAllDegeneratedTriangles(); 4649 4359 TriangleMap::iterator finder; 4650 4360 BoundaryTriangleSet *triangle = NULL, *partnerTriangle = NULL; 4651 int count = 0; 4652 4653 for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner) { 4361 int count = 0; 4362 4363 for (IndexToIndex::iterator TriangleKeyRunner = DegeneratedTriangles->begin(); 4364 TriangleKeyRunner != DegeneratedTriangles->end(); ++TriangleKeyRunner 4365 ) { 4654 4366 finder = TrianglesOnBoundary.find(TriangleKeyRunner->first); 4655 4367 if (finder != TrianglesOnBoundary.end()) … … 4668 4380 trianglesShareLine = trianglesShareLine || triangle->lines[i] == partnerTriangle->lines[j]; 4669 4381 4670 if (trianglesShareLine && (triangle->endpoints[1]->LinesCount > 2) && (triangle->endpoints[2]->LinesCount > 2) && (triangle->endpoints[0]->LinesCount > 2)) { 4382 if (trianglesShareLine 4383 && (triangle->endpoints[1]->LinesCount > 2) 4384 && (triangle->endpoints[2]->LinesCount > 2) 4385 && (triangle->endpoints[0]->LinesCount > 2) 4386 ) { 4671 4387 // check whether we have to fix lines 4672 4388 BoundaryTriangleSet *Othertriangle = NULL; … … 4688 4404 // the line of triangle receives the degenerated ones 4689 4405 triangle->lines[i]->triangles.erase(Othertriangle->Nr); 4690 triangle->lines[i]->triangles.insert( TrianglePair(partnerTriangle->Nr, partnerTriangle));4691 for (int k = 0; k < 3;k++)4406 triangle->lines[i]->triangles.insert( TrianglePair( partnerTriangle->Nr, partnerTriangle) ); 4407 for (int k=0;k<3;k++) 4692 4408 if (triangle->lines[i] == Othertriangle->lines[k]) { 4693 4409 Othertriangle->lines[k] = partnerTriangle->lines[j]; … … 4695 4411 } 4696 4412 // the line of partnerTriangle receives the non-degenerated ones 4697 partnerTriangle->lines[j]->triangles.erase( partnerTriangle->Nr);4698 partnerTriangle->lines[j]->triangles.insert( TrianglePair(Othertriangle->Nr, Othertriangle));4413 partnerTriangle->lines[j]->triangles.erase( partnerTriangle->Nr); 4414 partnerTriangle->lines[j]->triangles.insert( TrianglePair( Othertriangle->Nr, Othertriangle) ); 4699 4415 partnerTriangle->lines[j] = triangle->lines[i]; 4700 4416 } … … 4702 4418 // erase the pair 4703 4419 count += (int) DegeneratedTriangles->erase(triangle->Nr); 4704 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl);4420 Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *triangle << "." << endl; 4705 4421 RemoveTesselationTriangle(triangle); 4706 4422 count += (int) DegeneratedTriangles->erase(partnerTriangle->Nr); 4707 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl);4423 Log() << Verbose(0) << "RemoveDegeneratedTriangles() removes triangle " << *partnerTriangle << "." << endl; 4708 4424 RemoveTesselationTriangle(partnerTriangle); 4709 4425 } else { 4710 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle << " and its partner " << *partnerTriangle << " because it is essential for at" << " least one of the endpoints to be kept in the tesselation structure." << endl); 4711 } 4712 } 4713 delete (DegeneratedTriangles); 4426 Log() << Verbose(0) << "RemoveDegeneratedTriangles() does not remove triangle " << *triangle 4427 << " and its partner " << *partnerTriangle << " because it is essential for at" 4428 << " least one of the endpoints to be kept in the tesselation structure." << endl; 4429 } 4430 } 4431 delete(DegeneratedTriangles); 4714 4432 if (count > 0) 4715 4433 LastTriangle = NULL; 4716 4434 4717 DoLog(0) && (Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl);4435 Log() << Verbose(0) << "RemoveDegeneratedTriangles() removed " << count << " triangles:" << endl; 4718 4436 } 4719 4437 … … 4728 4446 void Tesselation::AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC) 4729 4447 { 4730 Info FunctionInfo(__func__);4448 Info FunctionInfo(__func__); 4731 4449 // find nearest boundary point 4732 4450 class TesselPoint *BackupPoint = NULL; … … 4741 4459 NearestBoundaryPoint = PointRunner->second; 4742 4460 } else { 4743 DoeLog(1) && (eLog() << Verbose(1) << "I cannot find the boundary point." << endl);4461 eLog() << Verbose(1) << "I cannot find the boundary point." << endl; 4744 4462 return; 4745 4463 } 4746 DoLog(0) && (Log() << Verbose(0) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl);4464 Log() << Verbose(0) << "Nearest point on boundary is " << NearestPoint->Name << "." << endl; 4747 4465 4748 4466 // go through its lines and find the best one to split … … 4759 4477 CenterToPoint.SubtractVector(point->node); 4760 4478 angle = CenterToPoint.Angle(&BaseLine); 4761 if (fabs(angle - M_PI / 2.) < fabs(BestAngle - M_PI /2.)) {4479 if (fabs(angle - M_PI/2.) < fabs(BestAngle - M_PI/2.)) { 4762 4480 BestAngle = angle; 4763 4481 BestLine = Runner->second; … … 4769 4487 BestLine->triangles.erase(TempTriangle->Nr); 4770 4488 int nr = -1; 4771 for (int i = 0; i <3; i++) {4489 for (int i=0;i<3; i++) { 4772 4490 if (TempTriangle->lines[i] == BestLine) { 4773 4491 nr = i; … … 4777 4495 4778 4496 // create new triangle to connect point (connects automatically with the missing spot of the chosen line) 4779 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);4497 Log() << Verbose(2) << "Adding new triangle points."<< endl; 4780 4498 AddTesselationPoint((BestLine->endpoints[0]->node), 0); 4781 4499 AddTesselationPoint((BestLine->endpoints[1]->node), 1); 4782 4500 AddTesselationPoint(point, 2); 4783 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);4784 AddTesselationLine( NULL, NULL,TPS[0], TPS[1], 0);4785 AddTesselationLine( NULL, NULL,TPS[0], TPS[2], 1);4786 AddTesselationLine( NULL, NULL,TPS[1], TPS[2], 2);4501 Log() << Verbose(2) << "Adding new triangle lines."<< endl; 4502 AddTesselationLine(TPS[0], TPS[1], 0); 4503 AddTesselationLine(TPS[0], TPS[2], 1); 4504 AddTesselationLine(TPS[1], TPS[2], 2); 4787 4505 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 4788 4506 BTS->GetNormalVector(TempTriangle->NormalVector); 4789 4507 BTS->NormalVector.Scale(-1.); 4790 DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl);4508 Log() << Verbose(1) << "INFO: NormalVector of new triangle is " << BTS->NormalVector << "." << endl; 4791 4509 AddTesselationTriangle(); 4792 4510 4793 4511 // create other side of this triangle and close both new sides of the first created triangle 4794 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle points." << endl);4512 Log() << Verbose(2) << "Adding new triangle points."<< endl; 4795 4513 AddTesselationPoint((BestLine->endpoints[0]->node), 0); 4796 4514 AddTesselationPoint((BestLine->endpoints[1]->node), 1); 4797 4515 AddTesselationPoint(point, 2); 4798 DoLog(2) && (Log() << Verbose(2) << "Adding new triangle lines." << endl);4799 AddTesselationLine( NULL, NULL,TPS[0], TPS[1], 0);4800 AddTesselationLine( NULL, NULL,TPS[0], TPS[2], 1);4801 AddTesselationLine( NULL, NULL,TPS[1], TPS[2], 2);4516 Log() << Verbose(2) << "Adding new triangle lines."<< endl; 4517 AddTesselationLine(TPS[0], TPS[1], 0); 4518 AddTesselationLine(TPS[0], TPS[2], 1); 4519 AddTesselationLine(TPS[1], TPS[2], 2); 4802 4520 BTS = new class BoundaryTriangleSet(BLS, TrianglesOnBoundaryCount); 4803 4521 BTS->GetNormalVector(TempTriangle->NormalVector); 4804 DoLog(1) && (Log() << Verbose(1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl);4522 Log() << Verbose(1) << "INFO: NormalVector of other new triangle is " << BTS->NormalVector << "." << endl; 4805 4523 AddTesselationTriangle(); 4806 4524 4807 4525 // add removed triangle to the last open line of the second triangle 4808 for (int i = 0; i < 3;i++) { // look for the same line as BestLine (only it's its degenerated companion)4526 for (int i=0;i<3;i++) { // look for the same line as BestLine (only it's its degenerated companion) 4809 4527 if ((BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[0])) && (BTS->lines[i]->ContainsBoundaryPoint(BestLine->endpoints[1]))) { 4810 if (BestLine == BTS->lines[i]) {4811 DoeLog(0) && (eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl);4528 if (BestLine == BTS->lines[i]){ 4529 eLog() << Verbose(0) << "BestLine is same as found line, something's wrong here!" << endl; 4812 4530 performCriticalExit(); 4813 4531 } 4814 BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle));4532 BTS->lines[i]->triangles.insert( pair<int, class BoundaryTriangleSet *> (TempTriangle->Nr, TempTriangle) ); 4815 4533 TempTriangle->lines[nr] = BTS->lines[i]; 4816 4534 break; 4817 4535 } 4818 4536 } 4819 } 4820 ; 4537 }; 4821 4538 4822 4539 /** Writes the envelope to file. … … 4827 4544 void Tesselation::Output(const char *filename, const PointCloud * const cloud) 4828 4545 { 4829 Info FunctionInfo(__func__);4546 Info FunctionInfo(__func__); 4830 4547 ofstream *tempstream = NULL; 4831 4548 string NameofTempFile; … … 4833 4550 4834 4551 if (LastTriangle != NULL) { 4835 sprintf(NumberName, "-%04d-%s_%s_%s", (int) TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name);4552 sprintf(NumberName, "-%04d-%s_%s_%s", (int)TrianglesOnBoundary.size(), LastTriangle->endpoints[0]->node->Name, LastTriangle->endpoints[1]->node->Name, LastTriangle->endpoints[2]->node->Name); 4836 4553 if (DoTecplotOutput) { 4837 4554 string NameofTempFile(filename); 4838 4555 NameofTempFile.append(NumberName); 4839 for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))4840 NameofTempFile.erase(npos, 1);4556 for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos)) 4557 NameofTempFile.erase(npos, 1); 4841 4558 NameofTempFile.append(TecplotSuffix); 4842 DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");4559 Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 4843 4560 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 4844 4561 WriteTecplotFile(tempstream, this, cloud, TriangleFilesWritten); 4845 4562 tempstream->close(); 4846 4563 tempstream->flush(); 4847 delete (tempstream);4564 delete(tempstream); 4848 4565 } 4849 4566 … … 4851 4568 string NameofTempFile(filename); 4852 4569 NameofTempFile.append(NumberName); 4853 for (size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos))4854 NameofTempFile.erase(npos, 1);4570 for(size_t npos = NameofTempFile.find_first_of(' '); npos != string::npos; npos = NameofTempFile.find(' ', npos)) 4571 NameofTempFile.erase(npos, 1); 4855 4572 NameofTempFile.append(Raster3DSuffix); 4856 DoLog(0) && (Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n");4573 Log() << Verbose(0) << "Writing temporary non convex hull to file " << NameofTempFile << ".\n"; 4857 4574 tempstream = new ofstream(NameofTempFile.c_str(), ios::trunc); 4858 4575 WriteRaster3dFile(tempstream, this, cloud); … … 4860 4577 tempstream->close(); 4861 4578 tempstream->flush(); 4862 delete (tempstream);4579 delete(tempstream); 4863 4580 } 4864 4581 } 4865 4582 if (DoTecplotOutput || DoRaster3DOutput) 4866 4583 TriangleFilesWritten++; 4867 } 4868 ; 4869 4870 struct BoundaryPolygonSetCompare 4871 { 4872 bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const 4873 { 4584 }; 4585 4586 struct BoundaryPolygonSetCompare { 4587 bool operator()(const BoundaryPolygonSet * s1, const BoundaryPolygonSet * s2) const { 4874 4588 if (s1->endpoints.size() < s2->endpoints.size()) 4875 4589 return true; … … 4900 4614 { 4901 4615 Info FunctionInfo(__func__); 4616 4902 4617 /// 2. Go through all BoundaryPointSet's, check their triangles' NormalVector 4903 4618 IndexToIndex *DegeneratedTriangles = FindAllDegeneratedTriangles(); 4904 set <BoundaryPointSet *> EndpointCandidateList;4905 pair <set<BoundaryPointSet *>::iterator, bool> InsertionTester;4906 pair <map<int, Vector *>::iterator, bool> TriangleInsertionTester;4619 set < BoundaryPointSet *> EndpointCandidateList; 4620 pair < set < BoundaryPointSet *>::iterator, bool > InsertionTester; 4621 pair < map < int, Vector *>::iterator, bool > TriangleInsertionTester; 4907 4622 for (PointMap::const_iterator Runner = PointsOnBoundary.begin(); Runner != PointsOnBoundary.end(); Runner++) { 4908 DoLog(0) && (Log() << Verbose(0) << "Current point is " << *Runner->second << "." << endl);4909 map <int, Vector *> TriangleVectors;4623 Log() << Verbose(0) << "Current point is " << *Runner->second << "." << endl; 4624 map < int, Vector *> TriangleVectors; 4910 4625 // gather all NormalVectors 4911 DoLog(1) && (Log() << Verbose(1) << "Gathering triangles ..." << endl);4626 Log() << Verbose(1) << "Gathering triangles ..." << endl; 4912 4627 for (LineMap::const_iterator LineRunner = (Runner->second)->lines.begin(); LineRunner != (Runner->second)->lines.end(); LineRunner++) 4913 4628 for (TriangleMap::const_iterator TriangleRunner = (LineRunner->second)->triangles.begin(); TriangleRunner != (LineRunner->second)->triangles.end(); TriangleRunner++) { 4914 4629 if (DegeneratedTriangles->find(TriangleRunner->second->Nr) == DegeneratedTriangles->end()) { 4915 TriangleInsertionTester = TriangleVectors.insert( pair<int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)));4630 TriangleInsertionTester = TriangleVectors.insert( pair< int, Vector *> ((TriangleRunner->second)->Nr, &((TriangleRunner->second)->NormalVector)) ); 4916 4631 if (TriangleInsertionTester.second) 4917 DoLog(1) && (Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl);4632 Log() << Verbose(1) << " Adding triangle " << *(TriangleRunner->second) << " to triangles to check-list." << endl; 4918 4633 } else { 4919 DoLog(1) && (Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl);4634 Log() << Verbose(1) << " NOT adding triangle " << *(TriangleRunner->second) << " as it's a simply degenerated one." << endl; 4920 4635 } 4921 4636 } 4922 4637 // check whether there are two that are parallel 4923 DoLog(1) && (Log() << Verbose(1) << "Finding two parallel triangles ..." << endl);4924 for (map <int, Vector *>::iterator VectorWalker = TriangleVectors.begin(); VectorWalker != TriangleVectors.end(); VectorWalker++)4925 for (map <int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++)4638 Log() << Verbose(1) << "Finding two parallel triangles ..." << endl; 4639 for (map < int, Vector *>::iterator VectorWalker = TriangleVectors.begin(); VectorWalker != TriangleVectors.end(); VectorWalker++) 4640 for (map < int, Vector *>::iterator VectorRunner = VectorWalker; VectorRunner != TriangleVectors.end(); VectorRunner++) 4926 4641 if (VectorWalker != VectorRunner) { // skip equals 4927 const double SCP = VectorWalker->second->ScalarProduct(VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles4928 DoLog(1) && (Log() << Verbose(1) << "Checking " << *VectorWalker->second << " against " << *VectorRunner->second << ": " << SCP << endl);4642 const double SCP = VectorWalker->second->ScalarProduct(VectorRunner->second); // ScalarProduct should result in -1. for degenerated triangles 4643 Log() << Verbose(1) << "Checking " << *VectorWalker->second<< " against " << *VectorRunner->second << ": " << SCP << endl; 4929 4644 if (fabs(SCP + 1.) < ParallelEpsilon) { 4930 4645 InsertionTester = EndpointCandidateList.insert((Runner->second)); 4931 4646 if (InsertionTester.second) 4932 DoLog(0) && (Log() << Verbose(0) << " Adding " << *Runner->second << " to endpoint candidate list." << endl);4647 Log() << Verbose(0) << " Adding " << *Runner->second << " to endpoint candidate list." << endl; 4933 4648 // and break out of both loops 4934 4649 VectorWalker = TriangleVectors.end(); … … 4938 4653 } 4939 4654 } 4940 delete (DegeneratedTriangles); 4655 4941 4656 /// 3. Find connected endpoint candidates and put them into a polygon 4942 4657 UniquePolygonSet ListofDegeneratedPolygons; … … 4944 4659 BoundaryPointSet *OtherWalker = NULL; 4945 4660 BoundaryPolygonSet *Current = NULL; 4946 stack <BoundaryPointSet*> ToCheckConnecteds;4661 stack <BoundaryPointSet*> ToCheckConnecteds; 4947 4662 while (!EndpointCandidateList.empty()) { 4948 4663 Walker = *(EndpointCandidateList.begin()); 4949 if (Current == NULL) { // create a new polygon with current candidate4950 DoLog(0) && (Log() << Verbose(0) << "Starting new polygon set at point " << *Walker << endl);4664 if (Current == NULL) { // create a new polygon with current candidate 4665 Log() << Verbose(0) << "Starting new polygon set at point " << *Walker << endl; 4951 4666 Current = new BoundaryPolygonSet; 4952 4667 Current->endpoints.insert(Walker); … … 4961 4676 for (LineMap::const_iterator LineWalker = Walker->lines.begin(); LineWalker != Walker->lines.end(); LineWalker++) { 4962 4677 OtherWalker = (LineWalker->second)->GetOtherEndpoint(Walker); 4963 DoLog(1) && (Log() << Verbose(1) << "Checking " << *OtherWalker << endl);4964 set <BoundaryPointSet *>::iterator Finder = EndpointCandidateList.find(OtherWalker);4965 if (Finder != EndpointCandidateList.end()) { // found a connected partner4966 DoLog(1) && (Log() << Verbose(1) << " Adding to polygon." << endl);4678 Log() << Verbose(1) << "Checking " << *OtherWalker << endl; 4679 set < BoundaryPointSet *>::iterator Finder = EndpointCandidateList.find(OtherWalker); 4680 if (Finder != EndpointCandidateList.end()) { // found a connected partner 4681 Log() << Verbose(1) << " Adding to polygon." << endl; 4967 4682 Current->endpoints.insert(OtherWalker); 4968 EndpointCandidateList.erase(Finder); // remove from candidates4969 ToCheckConnecteds.push(OtherWalker); // but check its partners too4683 EndpointCandidateList.erase(Finder); // remove from candidates 4684 ToCheckConnecteds.push(OtherWalker); // but check its partners too 4970 4685 } else { 4971 DoLog(1) && (Log() << Verbose(1) << " is not connected to " << *Walker << endl);4686 Log() << Verbose(1) << " is not connected to " << *Walker << endl; 4972 4687 } 4973 4688 } 4974 4689 } 4975 4690 4976 DoLog(0) && (Log() << Verbose(0) << "Final polygon is " << *Current << endl);4691 Log() << Verbose(0) << "Final polygon is " << *Current << endl; 4977 4692 ListofDegeneratedPolygons.insert(Current); 4978 4693 Current = NULL; … … 4981 4696 const int counter = ListofDegeneratedPolygons.size(); 4982 4697 4983 DoLog(0) && (Log() << Verbose(0) << "The following " << counter << " degenerated polygons have been found: " << endl);4698 Log() << Verbose(0) << "The following " << counter << " degenerated polygons have been found: " << endl; 4984 4699 for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++) 4985 DoLog(0) && (Log() << Verbose(0) << " " << **PolygonRunner << endl);4700 Log() << Verbose(0) << " " << **PolygonRunner << endl; 4986 4701 4987 4702 /// 4. Go through all these degenerated polygons 4988 4703 for (UniquePolygonSet::iterator PolygonRunner = ListofDegeneratedPolygons.begin(); PolygonRunner != ListofDegeneratedPolygons.end(); PolygonRunner++) { 4989 stack <int> TriangleNrs;4704 stack <int> TriangleNrs; 4990 4705 Vector NormalVector; 4991 4706 /// 4a. Gather all triangles of this polygon … … 4994 4709 // check whether number is bigger than 2, otherwise it's just a simply degenerated one and nothing to do. 4995 4710 if (T->size() == 2) { 4996 DoLog(1) && (Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl);4997 delete (T);4711 Log() << Verbose(1) << " Skipping degenerated polygon, is just a (already simply degenerated) triangle." << endl; 4712 delete(T); 4998 4713 continue; 4999 4714 } … … 5004 4719 // connections to either polygon ... 5005 4720 if (T->size() % 2 != 0) { 5006 DoeLog(0) && (eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl);4721 eLog() << Verbose(0) << " degenerated polygon contains an odd number of triangles, probably contains bridging non-degenerated ones, too!" << endl; 5007 4722 performCriticalExit(); 5008 4723 } 5009 TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator 4724 4725 TriangleSet::iterator TriangleWalker = T->begin(); // is the inner iterator 5010 4726 /// 4a. Get NormalVector for one side (this is "front") 5011 4727 NormalVector.CopyVector(&(*TriangleWalker)->NormalVector); 5012 DoLog(1) && (Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl);4728 Log() << Verbose(1) << "\"front\" defining triangle is " << **TriangleWalker << " and Normal vector of \"front\" side is " << NormalVector << endl; 5013 4729 TriangleWalker++; 5014 4730 TriangleSet::iterator TriangleSprinter = TriangleWalker; // is the inner advanced iterator … … 5019 4735 triangle = *TriangleWalker; 5020 4736 TriangleSprinter++; 5021 DoLog(1) && (Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl);4737 Log() << Verbose(1) << "Current triangle to test for removal: " << *triangle << endl; 5022 4738 if (triangle->NormalVector.ScalarProduct(&NormalVector) < 0) { // if from other side, then delete and remove from list 5023 DoLog(1) && (Log() << Verbose(1) << " Removing ... " << endl);4739 Log() << Verbose(1) << " Removing ... " << endl; 5024 4740 TriangleNrs.push(triangle->Nr); 5025 4741 T->erase(TriangleWalker); 5026 4742 RemoveTesselationTriangle(triangle); 5027 4743 } else 5028 DoLog(1) && (Log() << Verbose(1) << " Keeping ... " << endl);4744 Log() << Verbose(1) << " Keeping ... " << endl; 5029 4745 } 5030 4746 /// 4c. Copy all "front" triangles but with inverse NormalVector 5031 4747 TriangleWalker = T->begin(); 5032 while (TriangleWalker != T->end()) { // go through all front triangles5033 DoLog(1) && (Log() << Verbose(1) << " Re-creating triangle " << **TriangleWalker << " with NormalVector " << (*TriangleWalker)->NormalVector << endl);4748 while (TriangleWalker != T->end()) { // go through all front triangles 4749 Log() << Verbose(1) << " Re-creating triangle " << **TriangleWalker << " with NormalVector " << (*TriangleWalker)->NormalVector << endl; 5034 4750 for (int i = 0; i < 3; i++) 5035 4751 AddTesselationPoint((*TriangleWalker)->endpoints[i]->node, i); 5036 AddTesselationLine( NULL, NULL,TPS[0], TPS[1], 0);5037 AddTesselationLine( NULL, NULL,TPS[0], TPS[2], 1);5038 AddTesselationLine( NULL, NULL,TPS[1], TPS[2], 2);4752 AddTesselationLine(TPS[0], TPS[1], 0); 4753 AddTesselationLine(TPS[0], TPS[2], 1); 4754 AddTesselationLine(TPS[1], TPS[2], 2); 5039 4755 if (TriangleNrs.empty()) 5040 DoeLog(0) && (eLog() << Verbose(0) << "No more free triangle numbers!" << endl);4756 eLog() << Verbose(0) << "No more free triangle numbers!" << endl; 5041 4757 BTS = new BoundaryTriangleSet(BLS, TriangleNrs.top()); // copy triangle ... 5042 4758 AddTesselationTriangle(); // ... and add … … 5047 4763 } 5048 4764 if (!TriangleNrs.empty()) { 5049 DoeLog(0) && (eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl); 5050 } 5051 delete (T); // remove the triangleset 5052 } 4765 eLog() << Verbose(0) << "There have been less triangles created than removed!" << endl; 4766 } 4767 delete(T); // remove the triangleset 4768 } 4769 5053 4770 IndexToIndex * SimplyDegeneratedTriangles = FindAllDegeneratedTriangles(); 5054 DoLog(0) && (Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl);4771 Log() << Verbose(0) << "Final list of simply degenerated triangles found, containing " << SimplyDegeneratedTriangles->size() << " triangles:" << endl; 5055 4772 IndexToIndex::iterator it; 5056 4773 for (it = SimplyDegeneratedTriangles->begin(); it != SimplyDegeneratedTriangles->end(); it++) 5057 DoLog(0) && (Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl); 5058 delete (SimplyDegeneratedTriangles); 4774 Log() << Verbose(0) << (*it).first << " => " << (*it).second << endl; 4775 delete(SimplyDegeneratedTriangles); 4776 5059 4777 /// 5. exit 5060 4778 UniquePolygonSet::iterator PolygonRunner; 5061 4779 while (!ListofDegeneratedPolygons.empty()) { 5062 4780 PolygonRunner = ListofDegeneratedPolygons.begin(); 5063 delete (*PolygonRunner);4781 delete(*PolygonRunner); 5064 4782 ListofDegeneratedPolygons.erase(PolygonRunner); 5065 4783 } 5066 4784 5067 4785 return counter; 5068 } 5069 ; 4786 };
Note:
See TracChangeset
for help on using the changeset viewer.
