- File:
-
- 1 edited
-
molecuilder/src/tesselationhelpers.cpp (modified) (48 diffs)
Legend:
- Unmodified
- Added
- Removed
-
molecuilder/src/tesselationhelpers.cpp
rbe2997 r543ce4 8 8 #include <fstream> 9 9 10 #include "info.hpp"11 10 #include "linkedcell.hpp" 12 11 #include "log.hpp" … … 16 15 #include "verbose.hpp" 17 16 18 double DetGet(gsl_matrix * const A, const int inPlace) 19 { 20 Info FunctionInfo(__func__); 17 double DetGet(gsl_matrix * const A, const int inPlace) { 21 18 /* 22 19 inPlace = 1 => A is replaced with the LU decomposed copy. … … 48 45 void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS) 49 46 { 50 Info FunctionInfo(__func__);51 47 gsl_matrix *A = gsl_matrix_calloc(3,3); 52 48 double m11, m12, m13, m14; … … 81 77 82 78 if (fabs(m11) < MYEPSILON) 83 eLog() << Verbose( 1) << "three points are colinear." << endl;79 eLog() << Verbose(0) << "ERROR: three points are colinear." << endl; 84 80 85 81 center->x[0] = 0.5 * m12/ m11; … … 88 84 89 85 if (fabs(a.Distance(center) - RADIUS) > MYEPSILON) 90 eLog() << Verbose( 1) << "The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl;86 eLog() << Verbose(0) << "ERROR: The given center is further way by " << fabs(a.Distance(center) - RADIUS) << " from a than RADIUS." << endl; 91 87 92 88 gsl_matrix_free(A); … … 115 111 const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius) 116 112 { 117 Info FunctionInfo(__func__);118 113 Vector TempNormal, helper; 119 114 double Restradius; 120 115 Vector OtherCenter; 116 Log() << Verbose(3) << "Begin of GetCenterOfSphere.\n"; 121 117 Center->Zero(); 122 118 helper.CopyVector(&a); … … 132 128 Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma))); 133 129 NewUmkreismittelpunkt->CopyVector(Center); 134 Log() << Verbose( 1) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n";130 Log() << Verbose(4) << "Center of new circumference is " << *NewUmkreismittelpunkt << ".\n"; 135 131 // Here we calculated center of circumscribing circle, using barycentric coordinates 136 Log() << Verbose( 1) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n";132 Log() << Verbose(4) << "Center of circumference is " << *Center << " in direction " << *Direction << ".\n"; 137 133 138 134 TempNormal.CopyVector(&a); … … 158 154 TempNormal.Normalize(); 159 155 Restradius = sqrt(RADIUS*RADIUS - Umkreisradius*Umkreisradius); 160 Log() << Verbose( 1) << "Height of center of circumference to center of sphere is " << Restradius << ".\n";156 Log() << Verbose(4) << "Height of center of circumference to center of sphere is " << Restradius << ".\n"; 161 157 TempNormal.Scale(Restradius); 162 Log() << Verbose( 1) << "Shift vector to sphere of circumference is " << TempNormal << ".\n";158 Log() << Verbose(4) << "Shift vector to sphere of circumference is " << TempNormal << ".\n"; 163 159 164 160 Center->AddVector(&TempNormal); 165 Log() << Verbose( 1) << "Center of sphere of circumference is " << *Center << ".\n";161 Log() << Verbose(0) << "Center of sphere of circumference is " << *Center << ".\n"; 166 162 GetSphere(&OtherCenter, a, b, c, RADIUS); 167 Log() << Verbose(1) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n"; 163 Log() << Verbose(0) << "OtherCenter of sphere of circumference is " << OtherCenter << ".\n"; 164 Log() << Verbose(3) << "End of GetCenterOfSphere.\n"; 168 165 }; 169 166 … … 177 174 void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c) 178 175 { 179 Info FunctionInfo(__func__);180 176 Vector helper; 181 177 double alpha, beta, gamma; … … 190 186 beta = M_PI - SideC.Angle(&SideA); 191 187 gamma = M_PI - SideA.Angle(&SideB); 192 //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl; 193 if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) { 194 eLog() << Verbose(1) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl; 195 } 188 //Log() << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl; 189 if (fabs(M_PI - alpha - beta - gamma) > HULLEPSILON) 190 eLog() << Verbose(0) << "GetCenterofCircumcircle: Sum of angles " << (alpha+beta+gamma)/M_PI*180. << " > 180 degrees by " << fabs(M_PI - alpha - beta - gamma)/M_PI*180. << "!" << endl; 196 191 197 192 Center->Zero(); … … 223 218 double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection) 224 219 { 225 Info FunctionInfo(__func__);226 220 Vector helper; 227 221 double radius, alpha; … … 230 224 // test whether new center is on the parameter circle's plane 231 225 if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) { 232 eLog() << Verbose( 1) << "Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal)) << "!" << endl;226 eLog() << Verbose(0) << "ERROR: Something's very wrong here: NewSphereCenter is not on the band's plane as desired by " <<fabs(helper.ScalarProduct(&CirclePlaneNormal)) << "!" << endl; 233 227 helper.ProjectOntoPlane(&CirclePlaneNormal); 234 228 } … … 236 230 // test whether the new center vector has length of CircleRadius 237 231 if (fabs(radius - CircleRadius) > HULLEPSILON) 238 eLog() << Verbose(1) << " The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;232 eLog() << Verbose(1) << "ERROR: The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl; 239 233 alpha = helper.Angle(&OldSphereCenter); 240 234 // make the angle unique by checking the halfplanes/search direction 241 235 if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON) // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals 242 236 alpha = 2.*M_PI - alpha; 243 //Log() << Verbose( 1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;237 //Log() << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl; 244 238 radius = helper.Distance(&OldSphereCenter); 245 239 helper.ProjectOntoPlane(&NormalVector); 246 240 // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles 247 241 if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) { 248 //Log() << Verbose( 1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;242 //Log() << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl; 249 243 return alpha; 250 244 } else { … … 269 263 double MinIntersectDistance(const gsl_vector * x, void *params) 270 264 { 271 Info FunctionInfo(__func__);272 265 double retval = 0; 273 266 struct Intersection *I = (struct Intersection *)params; … … 290 283 291 284 retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB); 292 //Log() << Verbose( 1) << "MinIntersectDistance called, result: " << retval << endl;285 //Log() << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl; 293 286 294 287 return retval; … … 310 303 bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4) 311 304 { 312 Info FunctionInfo(__func__);313 305 bool result; 314 306 … … 358 350 359 351 if (status == GSL_SUCCESS) { 360 Log() << Verbose( 1) << "converged to minimum" << endl;352 Log() << Verbose(2) << "converged to minimum" << endl; 361 353 } 362 354 } while (status == GSL_CONTINUE && iter < 100); … … 383 375 t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB); 384 376 385 Log() << Verbose( 1) << "Intersection " << intersection << " is at "377 Log() << Verbose(2) << "Intersection " << intersection << " is at " 386 378 << t1 << " for (" << point1 << "," << point2 << ") and at " 387 379 << t2 << " for (" << point3 << "," << point4 << "): "; 388 380 389 381 if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) { 390 Log() << Verbose( 1) << "true intersection." << endl;382 Log() << Verbose(0) << "true intersection." << endl; 391 383 result = true; 392 384 } else { 393 Log() << Verbose( 1) << "intersection out of region of interest." << endl;385 Log() << Verbose(0) << "intersection out of region of interest." << endl; 394 386 result = false; 395 387 } … … 414 406 double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector) 415 407 { 416 Info FunctionInfo(__func__);417 408 if (reference.IsZero()) 418 409 return M_PI; … … 426 417 } 427 418 428 Log() << Verbose( 1) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;419 Log() << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl; 429 420 430 421 return phi; … … 441 432 double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d) 442 433 { 443 Info FunctionInfo(__func__);444 434 Vector Point, TetraederVector[3]; 445 435 double volume; … … 465 455 bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3]) 466 456 { 467 Info FunctionInfo(__func__);468 457 bool result = false; 469 458 int counter = 0; … … 472 461 for (int i=0;i<3;i++) 473 462 for (int j=i+1; j<3; j++) { 474 if (nodes[i] == NULL) { 475 Log() << Verbose(1) << "Node nr. " << i << " is not yet present." << endl; 476 result = true; 477 } else if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) { // there already is a line 463 if (nodes[i]->lines.find(nodes[j]->node->nr) != nodes[i]->lines.end()) { // there already is a line 478 464 LineMap::const_iterator FindLine; 479 465 pair<LineMap::const_iterator,LineMap::const_iterator> FindPair; … … 492 478 } 493 479 if ((!result) && (counter > 1)) { 494 Log() << Verbose( 1) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;480 Log() << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl; 495 481 result = true; 496 482 } … … 499 485 500 486 501 ///** Sort function for the candidate list. 502 // */ 503 //bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2) 504 //{ 505 // Info FunctionInfo(__func__); 506 // Vector BaseLineVector, OrthogonalVector, helper; 507 // if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check 508 // eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl; 509 // //return false; 510 // exit(1); 511 // } 512 // // create baseline vector 513 // BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node); 514 // BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 515 // BaseLineVector.Normalize(); 516 // 517 // // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!) 518 // helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node); 519 // helper.SubtractVector(candidate1->point->node); 520 // OrthogonalVector.CopyVector(&helper); 521 // helper.VectorProduct(&BaseLineVector); 522 // OrthogonalVector.SubtractVector(&helper); 523 // OrthogonalVector.Normalize(); 524 // 525 // // calculate both angles and correct with in-plane vector 526 // helper.CopyVector(candidate1->point->node); 527 // helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 528 // double phi = BaseLineVector.Angle(&helper); 529 // if (OrthogonalVector.ScalarProduct(&helper) > 0) { 530 // phi = 2.*M_PI - phi; 531 // } 532 // helper.CopyVector(candidate2->point->node); 533 // helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 534 // double psi = BaseLineVector.Angle(&helper); 535 // if (OrthogonalVector.ScalarProduct(&helper) > 0) { 536 // psi = 2.*M_PI - psi; 537 // } 538 // 539 // Log() << Verbose(1) << *candidate1->point << " has angle " << phi << endl; 540 // Log() << Verbose(1) << *candidate2->point << " has angle " << psi << endl; 541 // 542 // // return comparison 543 // return phi < psi; 544 //}; 487 /** Sort function for the candidate list. 488 */ 489 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2) 490 { 491 Vector BaseLineVector, OrthogonalVector, helper; 492 if (candidate1->BaseLine != candidate2->BaseLine) { // sanity check 493 Log() << Verbose(0) << "ERROR: sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl; 494 //return false; 495 exit(1); 496 } 497 // create baseline vector 498 BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node); 499 BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 500 BaseLineVector.Normalize(); 501 502 // create normal in-plane vector to cope with acos() non-uniqueness on [0,2pi] (note that is pointing in the "right" direction already, hence ">0" test!) 503 helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node); 504 helper.SubtractVector(candidate1->point->node); 505 OrthogonalVector.CopyVector(&helper); 506 helper.VectorProduct(&BaseLineVector); 507 OrthogonalVector.SubtractVector(&helper); 508 OrthogonalVector.Normalize(); 509 510 // calculate both angles and correct with in-plane vector 511 helper.CopyVector(candidate1->point->node); 512 helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 513 double phi = BaseLineVector.Angle(&helper); 514 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 515 phi = 2.*M_PI - phi; 516 } 517 helper.CopyVector(candidate2->point->node); 518 helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node); 519 double psi = BaseLineVector.Angle(&helper); 520 if (OrthogonalVector.ScalarProduct(&helper) > 0) { 521 psi = 2.*M_PI - psi; 522 } 523 524 Log() << Verbose(2) << *candidate1->point << " has angle " << phi << endl; 525 Log() << Verbose(2) << *candidate2->point << " has angle " << psi << endl; 526 527 // return comparison 528 return phi < psi; 529 }; 545 530 546 531 /** … … 554 539 TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC) 555 540 { 556 Info FunctionInfo(__func__);557 541 TesselPoint* closestPoint = NULL; 558 542 TesselPoint* secondClosestPoint = NULL; … … 565 549 for(int i=0;i<NDIM;i++) // store indices of this cell 566 550 N[i] = LC->n[i]; 567 Log() << Verbose( 1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;551 Log() << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 568 552 569 553 LC->GetNeighbourBounds(Nlower, Nupper); 570 //Log() << Verbose( 1) << endl;554 //Log() << Verbose(0) << endl; 571 555 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 572 556 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 573 557 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 574 558 const LinkedNodes *List = LC->GetCurrentCell(); 575 //Log() << Verbose( 1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;559 //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl; 576 560 if (List != NULL) { 577 561 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 590 574 } 591 575 } else { 592 eLog() << Verbose( 1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","576 eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," 593 577 << LC->n[2] << " is invalid!" << endl; 594 578 } … … 609 593 TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC) 610 594 { 611 Info FunctionInfo(__func__);612 595 TesselPoint* closestPoint = NULL; 613 596 SecondPoint = NULL; … … 620 603 for(int i=0;i<NDIM;i++) // store indices of this cell 621 604 N[i] = LC->n[i]; 622 Log() << Verbose( 1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;605 Log() << Verbose(3) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl; 623 606 624 607 LC->GetNeighbourBounds(Nlower, Nupper); 625 //Log() << Verbose( 1) << endl;608 //Log() << Verbose(0) << endl; 626 609 for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++) 627 610 for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++) 628 611 for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) { 629 612 const LinkedNodes *List = LC->GetCurrentCell(); 630 //Log() << Verbose( 1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;613 //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl; 631 614 if (List != NULL) { 632 615 for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) { … … 639 622 distance = currentNorm; 640 623 closestPoint = (*Runner); 641 //Log() << Verbose( 1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;624 //Log() << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl; 642 625 } else if (currentNorm < secondDistance) { 643 626 secondDistance = currentNorm; 644 627 SecondPoint = (*Runner); 645 //Log() << Verbose( 1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;628 //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl; 646 629 } 647 630 } 648 631 } else { 649 eLog() << Verbose( 1) << "The current cell " << LC->n[0] << "," << LC->n[1] << ","632 eLog() << Verbose(0) << "ERROR: The current cell " << LC->n[0] << "," << LC->n[1] << "," 650 633 << LC->n[2] << " is invalid!" << endl; 651 634 } … … 653 636 // output 654 637 if (closestPoint != NULL) { 655 Log() << Verbose( 1) << "Closest point is " << *closestPoint;638 Log() << Verbose(2) << "Closest point is " << *closestPoint; 656 639 if (SecondPoint != NULL) 657 640 Log() << Verbose(0) << " and second closest is " << *SecondPoint; … … 669 652 Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase) 670 653 { 671 Info FunctionInfo(__func__);672 654 // construct the plane of the two baselines (i.e. take both their directional vectors) 673 655 Vector Normal; … … 680 662 Normal.VectorProduct(&OtherBaseline); 681 663 Normal.Normalize(); 682 Log() << Verbose( 1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;664 Log() << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl; 683 665 684 666 // project one offset point of OtherBase onto this plane (and add plane offset vector) … … 697 679 Normal.CopyVector(Intersection); 698 680 Normal.SubtractVector(Base->endpoints[0]->node->node); 699 Log() << Verbose( 1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;681 Log() << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl; 700 682 701 683 return Intersection; … … 710 692 double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle) 711 693 { 712 Info FunctionInfo(__func__);713 694 double distance = 0.; 714 695 if (x == NULL) { … … 727 708 void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud) 728 709 { 729 Info FunctionInfo(__func__);730 710 TesselPoint *Walker = NULL; 731 711 int i; … … 758 738 } 759 739 } else { 760 eLog() << Verbose( 1) << "Given vrmlfile is " << vrmlfile << "." << endl;740 eLog() << Verbose(0) << "ERROR: Given vrmlfile is " << vrmlfile << "." << endl; 761 741 } 762 742 delete(center); … … 771 751 void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud) 772 752 { 773 Info FunctionInfo(__func__);774 753 Vector helper; 775 776 if (Tess->LastTriangle != NULL) { 777 // include the current position of the virtual sphere in the temporary raster3d file 778 Vector *center = cloud->GetCenter(); 779 // make the circumsphere's center absolute again 780 helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node); 781 helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node); 782 helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node); 783 helper.Scale(1./3.); 784 helper.SubtractVector(center); 785 // and add to file plus translucency object 786 *rasterfile << "# current virtual sphere\n"; 787 *rasterfile << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n"; 788 *rasterfile << "2\n " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n"; 789 *rasterfile << "9\n terminating special property\n"; 790 delete(center); 791 } 754 // include the current position of the virtual sphere in the temporary raster3d file 755 Vector *center = cloud->GetCenter(); 756 // make the circumsphere's center absolute again 757 helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node); 758 helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node); 759 helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node); 760 helper.Scale(1./3.); 761 helper.SubtractVector(center); 762 // and add to file plus translucency object 763 *rasterfile << "# current virtual sphere\n"; 764 *rasterfile << "8\n 25.0 0.6 -1.0 -1.0 -1.0 0.2 0 0 0 0\n"; 765 *rasterfile << "2\n " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n"; 766 *rasterfile << "9\n terminating special property\n"; 767 delete(center); 792 768 }; 793 769 … … 800 776 void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud) 801 777 { 802 Info FunctionInfo(__func__);803 778 TesselPoint *Walker = NULL; 804 779 int i; … … 833 808 *rasterfile << "9\n# terminating special property\n"; 834 809 } else { 835 eLog() << Verbose( 1) << "Given rasterfile is " << rasterfile << "." << endl;810 eLog() << Verbose(0) << "ERROR: Given rasterfile is " << rasterfile << "." << endl; 836 811 } 837 812 IncludeSphereinRaster3D(rasterfile, Tess, cloud); … … 846 821 void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N) 847 822 { 848 Info FunctionInfo(__func__);849 823 if ((tecplot != NULL) && (TesselStruct != NULL)) { 850 824 // write header 851 825 *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl; 852 826 *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl; 853 *tecplot << "ZONE T=\""; 854 if (N < 0) { 855 *tecplot << cloud->GetName(); 856 } else { 857 *tecplot << N << "-"; 858 for (int i=0;i<3;i++) 859 *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name; 860 } 827 *tecplot << "ZONE T=\"" << N << "-"; 828 for (int i=0;i<3;i++) 829 *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name; 861 830 *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl; 862 831 int i=0; … … 867 836 868 837 // print atom coordinates 838 Log() << Verbose(2) << "The following triangles were created:"; 869 839 int Counter = 1; 870 840 TesselPoint *Walker = NULL; … … 876 846 *tecplot << endl; 877 847 // print connectivity 878 Log() << Verbose(1) << "The following triangles were created:" << endl;879 848 for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) { 880 Log() << Verbose( 1) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name << endl;849 Log() << Verbose(0) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name; 881 850 *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl; 882 851 } 883 852 delete[] (LookupList); 853 Log() << Verbose(0) << endl; 884 854 } 885 855 }; … … 892 862 void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct) 893 863 { 894 Info FunctionInfo(__func__);895 864 class BoundaryPointSet *point = NULL; 896 865 class BoundaryLineSet *line = NULL; 897 866 867 //Log() << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl; 898 868 // calculate remaining concavity 899 869 for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) { … … 903 873 for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) { 904 874 line = LineRunner->second; 905 //Log() << Verbose( 1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;875 //Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl; 906 876 if (!line->CheckConvexityCriterion()) 907 877 point->value += 1; 908 878 } 909 879 } 880 //Log() << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl; 910 881 }; 911 882 … … 918 889 bool CheckListOfBaselines(const Tesselation * const TesselStruct) 919 890 { 920 Info FunctionInfo(__func__);921 891 LineMap::const_iterator testline; 922 892 bool result = false; … … 926 896 for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) { 927 897 if (testline->second->triangles.size() != 2) { 928 Log() << Verbose( 2) << *testline->second << "\t" << testline->second->triangles.size() << endl;898 Log() << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl; 929 899 counter++; 930 900 }
Note:
See TracChangeset
for help on using the changeset viewer.
