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