Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/tesselationhelpers.cpp

    re359a8 rb998c3  
    88#include <fstream>
    99
     10#include "info.hpp"
    1011#include "linkedcell.hpp"
    1112#include "log.hpp"
    … …  
    1516#include "verbose.hpp"
    1617
    17 double DetGet(gsl_matrix * const A, const int inPlace) {
     18double DetGet(gsl_matrix * const A, const int inPlace)
     19{
     20        Info FunctionInfo(__func__);
    1821  /*
    1922  inPlace = 1 => A is replaced with the LU decomposed copy.
    … …  
    4548void GetSphere(Vector * const center, const Vector &a, const Vector &b, const Vector &c, const double RADIUS)
    4649{
     50        Info FunctionInfo(__func__);
    4751  gsl_matrix *A = gsl_matrix_calloc(3,3);
    4852  double m11, m12, m13, m14;
    … …  
    111115    const double HalfplaneIndicator, const double AlternativeIndicator, const double alpha, const double beta, const double gamma, const double RADIUS, const double Umkreisradius)
    112116{
     117        Info FunctionInfo(__func__);
    113118  Vector TempNormal, helper;
    114119  double Restradius;
    115120  Vector OtherCenter;
    116   Log() << Verbose(3) << "Begin of GetCenterOfSphere.\n";
    117121  Center->Zero();
    118122  helper.CopyVector(&a);
    … …  
    128132  Center->Scale(1./(sin(2.*alpha) + sin(2.*beta) + sin(2.*gamma)));
    129133  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";
    131135  // 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";
    133137
    134138  TempNormal.CopyVector(&a);
    … …  
    154158  TempNormal.Normalize();
    155159  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";
    157161  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";
    159163
    160164  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";
    162166  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";
    165168};
    166169
    … …  
    174177void GetCenterofCircumcircle(Vector * const Center, const Vector &a, const Vector &b, const Vector &c)
    175178{
     179        Info FunctionInfo(__func__);
    176180  Vector helper;
    177181  double alpha, beta, gamma;
    … …  
    186190  beta = M_PI - SideC.Angle(&SideA);
    187191  gamma = M_PI - SideA.Angle(&SideB);
    188   //Log() << Verbose(3) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
     192  //Log() << Verbose(1) << "INFO: alpha = " << alpha/M_PI*180. << ", beta = " << beta/M_PI*180. << ", gamma = " << gamma/M_PI*180. << "." << endl;
    189193  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;
    191     performCriticalExit();
     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;
    192195  }
    193196
    … …  
    220223double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection)
    221224{
     225        Info FunctionInfo(__func__);
    222226  Vector helper;
    223227  double radius, alpha;
    224 
    225   helper.CopyVector(&NewSphereCenter);
     228  Vector RelativeOldSphereCenter;
     229  Vector RelativeNewSphereCenter;
     230
     231  RelativeOldSphereCenter.CopyVector(&OldSphereCenter);
     232  RelativeOldSphereCenter.SubtractVector(&CircleCenter);
     233  RelativeNewSphereCenter.CopyVector(&NewSphereCenter);
     234  RelativeNewSphereCenter.SubtractVector(&CircleCenter);
     235  helper.CopyVector(&RelativeNewSphereCenter);
    226236  // test whether new center is on the parameter circle's plane
    227237  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    … …  
    229239    helper.ProjectOntoPlane(&CirclePlaneNormal);
    230240  }
    231   radius = helper.ScalarProduct(&helper);
     241  radius = helper.NormSquared();
    232242  // test whether the new center vector has length of CircleRadius
    233243  if (fabs(radius - CircleRadius) > HULLEPSILON)
    234244    eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    235   alpha = helper.Angle(&OldSphereCenter);
     245  alpha = helper.Angle(&RelativeOldSphereCenter);
    236246  // make the angle unique by checking the halfplanes/search direction
    237247  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    238248    alpha = 2.*M_PI - alpha;
    239   //Log() << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
    240   radius = helper.Distance(&OldSphereCenter);
     249  Log() << Verbose(1) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << RelativeOldSphereCenter << " and resulting angle is " << alpha << "." << endl;
     250  radius = helper.Distance(&RelativeOldSphereCenter);
    241251  helper.ProjectOntoPlane(&NormalVector);
    242252  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
    243253  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
    244     //Log() << Verbose(2) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
     254    Log() << Verbose(1) << "INFO: Distance between old and new center is " << radius << " and between new center and baseline center is " << helper.Norm() << "." << endl;
    245255    return alpha;
    246256  } else {
    247     //Log() << Verbose(1) << "INFO: NewSphereCenter " << helper << " is too close to OldSphereCenter" << OldSphereCenter << "." << endl;
     257    Log() << Verbose(1) << "INFO: NewSphereCenter " << RelativeNewSphereCenter << " is too close to RelativeOldSphereCenter" << RelativeOldSphereCenter << "." << endl;
    248258    return 2.*M_PI;
    249259  }
    … …  
    265275double MinIntersectDistance(const gsl_vector * x, void *params)
    266276{
     277        Info FunctionInfo(__func__);
    267278  double retval = 0;
    268279  struct Intersection *I = (struct Intersection *)params;
    … …  
    285296
    286297  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
    287   //Log() << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
     298  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
    288299
    289300  return retval;
    … …  
    305316bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
    306317{
     318        Info FunctionInfo(__func__);
    307319  bool result;
    308320
    … …  
    352364
    353365        if (status == GSL_SUCCESS) {
    354           Log() << Verbose(2) << "converged to minimum" <<  endl;
     366          Log() << Verbose(1) << "converged to minimum" <<  endl;
    355367        }
    356368    } while (status == GSL_CONTINUE && iter < 100);
    … …  
    377389  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
    378390
    379   Log() << Verbose(2) << "Intersection " << intersection << " is at "
     391  Log() << Verbose(1) << "Intersection " << intersection << " is at "
    380392    << t1 << " for (" << point1 << "," << point2 << ") and at "
    381393    << t2 << " for (" << point3 << "," << point4 << "): ";
    382394
    383395  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
    384     Log() << Verbose(0) << "true intersection." << endl;
     396    Log() << Verbose(1) << "true intersection." << endl;
    385397    result = true;
    386398  } else {
    387     Log() << Verbose(0) << "intersection out of region of interest." << endl;
     399    Log() << Verbose(1) << "intersection out of region of interest." << endl;
    388400    result = false;
    389401  }
    … …  
    408420double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector)
    409421{
     422        Info FunctionInfo(__func__);
    410423  if (reference.IsZero())
    411424    return M_PI;
    … …  
    419432  }
    420433
    421   Log() << Verbose(4) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
     434  Log() << Verbose(1) << "INFO: " << point << " has angle " << phi << " with respect to reference " << reference << "." << endl;
    422435
    423436  return phi;
    … …  
    434447double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
    435448{
     449        Info FunctionInfo(__func__);
    436450  Vector Point, TetraederVector[3];
    437451  double volume;
    … …  
    457471bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3])
    458472{
     473        Info FunctionInfo(__func__);
    459474  bool result = false;
    460475  int counter = 0;
    … …  
    483498    }
    484499  if ((!result) && (counter > 1)) {
    485     Log() << Verbose(2) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
     500    Log() << Verbose(1) << "INFO: Degenerate triangle is ok, at least two, here " << counter << ", existing lines are used." << endl;
    486501    result = true;
    487502  }
    … …  
    490505
    491506
    492 /** Sort function for the candidate list.
    493  */
    494 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
    495 {
    496   Vector BaseLineVector, OrthogonalVector, helper;
    497   if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
    498     eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
    499     //return false;
    500     exit(1);
    501   }
    502   // create baseline vector
    503   BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
    504   BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    505   BaseLineVector.Normalize();
    506 
    507   // 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!)
    508   helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
    509   helper.SubtractVector(candidate1->point->node);
    510   OrthogonalVector.CopyVector(&helper);
    511   helper.VectorProduct(&BaseLineVector);
    512   OrthogonalVector.SubtractVector(&helper);
    513   OrthogonalVector.Normalize();
    514 
    515   // calculate both angles and correct with in-plane vector
    516   helper.CopyVector(candidate1->point->node);
    517   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    518   double phi = BaseLineVector.Angle(&helper);
    519   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    520     phi = 2.*M_PI - phi;
    521   }
    522   helper.CopyVector(candidate2->point->node);
    523   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    524   double psi = BaseLineVector.Angle(&helper);
    525   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    526     psi = 2.*M_PI - psi;
    527   }
    528 
    529   Log() << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
    530   Log() << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
    531 
    532   // return comparison
    533   return phi < psi;
    534 };
     507///** Sort function for the candidate list.
     508// */
     509//bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
     510//{
     511//      Info FunctionInfo(__func__);
     512//  Vector BaseLineVector, OrthogonalVector, helper;
     513//  if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
     514//    eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
     515//    //return false;
     516//    exit(1);
     517//  }
     518//  // create baseline vector
     519//  BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
     520//  BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     521//  BaseLineVector.Normalize();
     522//
     523//  // 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!)
     524//  helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
     525//  helper.SubtractVector(candidate1->point->node);
     526//  OrthogonalVector.CopyVector(&helper);
     527//  helper.VectorProduct(&BaseLineVector);
     528//  OrthogonalVector.SubtractVector(&helper);
     529//  OrthogonalVector.Normalize();
     530//
     531//  // calculate both angles and correct with in-plane vector
     532//  helper.CopyVector(candidate1->point->node);
     533//  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     534//  double phi = BaseLineVector.Angle(&helper);
     535//  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
     536//    phi = 2.*M_PI - phi;
     537//  }
     538//  helper.CopyVector(candidate2->point->node);
     539//  helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
     540//  double psi = BaseLineVector.Angle(&helper);
     541//  if (OrthogonalVector.ScalarProduct(&helper) > 0) {
     542//    psi = 2.*M_PI - psi;
     543//  }
     544//
     545//  Log() << Verbose(1) << *candidate1->point << " has angle " << phi << endl;
     546//  Log() << Verbose(1) << *candidate2->point << " has angle " << psi << endl;
     547//
     548//  // return comparison
     549//  return phi < psi;
     550//};
    535551
    536552/**
    … …  
    544560TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC)
    545561{
     562        Info FunctionInfo(__func__);
    546563  TesselPoint* closestPoint = NULL;
    547564  TesselPoint* secondClosestPoint = NULL;
    … …  
    554571  for(int i=0;i<NDIM;i++) // store indices of this cell
    555572    N[i] = LC->n[i];
    556   Log() << Verbose(2) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     573  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    557574
    558575  LC->GetNeighbourBounds(Nlower, Nupper);
    559   //Log() << Verbose(0) << endl;
     576  //Log() << Verbose(1) << endl;
    560577  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    561578    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    562579      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    563580        const LinkedNodes *List = LC->GetCurrentCell();
    564         //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     581        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    565582        if (List != NULL) {
    566583          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    … …  
    598615TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
    599616{
     617        Info FunctionInfo(__func__);
    600618  TesselPoint* closestPoint = NULL;
    601619  SecondPoint = NULL;
    … …  
    608626  for(int i=0;i<NDIM;i++) // store indices of this cell
    609627    N[i] = LC->n[i];
    610   Log() << Verbose(3) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
     628  Log() << Verbose(1) << "INFO: Center cell is " << N[0] << ", " << N[1] << ", " << N[2] << " with No. " << LC->index << "." << endl;
    611629
    612630  LC->GetNeighbourBounds(Nlower, Nupper);
    613   //Log() << Verbose(0) << endl;
     631  //Log() << Verbose(1) << endl;
    614632  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    615633    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    616634      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    617635        const LinkedNodes *List = LC->GetCurrentCell();
    618         //Log() << Verbose(3) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
     636        //Log() << Verbose(1) << "The current cell " << LC->n[0] << "," << LC->n[1] << "," << LC->n[2] << endl;
    619637        if (List != NULL) {
    620638          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    … …  
    627645              distance = currentNorm;
    628646              closestPoint = (*Runner);
    629               //Log() << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
     647              //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
    630648            } else if (currentNorm < secondDistance) {
    631649              secondDistance = currentNorm;
    632650              SecondPoint = (*Runner);
    633               //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
     651              //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
    634652            }
    635653          }
    … …  
    641659  // output
    642660  if (closestPoint != NULL) {
    643     Log() << Verbose(2) << "Closest point is " << *closestPoint;
     661    Log() << Verbose(1) << "Closest point is " << *closestPoint;
    644662    if (SecondPoint != NULL)
    645663      Log() << Verbose(0) << " and second closest is " << *SecondPoint;
    … …  
    657675Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase)
    658676{
     677        Info FunctionInfo(__func__);
    659678  // construct the plane of the two baselines (i.e. take both their directional vectors)
    660679  Vector Normal;
    … …  
    667686  Normal.VectorProduct(&OtherBaseline);
    668687  Normal.Normalize();
    669   Log() << Verbose(4) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
     688  Log() << Verbose(1) << "First direction is " << Baseline << ", second direction is " << OtherBaseline << ", normal of intersection plane is " << Normal << "." << endl;
    670689
    671690  // project one offset point of OtherBase onto this plane (and add plane offset vector)
    … …  
    684703  Normal.CopyVector(Intersection);
    685704  Normal.SubtractVector(Base->endpoints[0]->node->node);
    686   Log() << Verbose(3) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
     705  Log() << Verbose(1) << "Found closest point on " << *Base << " at " << *Intersection << ", factor in line is " << fabs(Normal.ScalarProduct(&Baseline)/Baseline.NormSquared()) << "." << endl;
    687706
    688707  return Intersection;
    … …  
    697716double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle)
    698717{
     718        Info FunctionInfo(__func__);
    699719  double distance = 0.;
    700720  if (x == NULL) {
    … …  
    713733void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud)
    714734{
     735        Info FunctionInfo(__func__);
    715736  TesselPoint *Walker = NULL;
    716737  int i;
    … …  
    756777void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    757778{
     779        Info FunctionInfo(__func__);
    758780  Vector helper;
    759   // include the current position of the virtual sphere in the temporary raster3d file
    760   Vector *center = cloud->GetCenter();
    761   // make the circumsphere's center absolute again
    762   helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
    763   helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
    764   helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
    765   helper.Scale(1./3.);
    766   helper.SubtractVector(center);
    767   // and add to file plus translucency object
    768   *rasterfile << "# current virtual sphere\n";
    769   *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
    770   *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
    771   *rasterfile << "9\n  terminating special property\n";
    772   delete(center);
     781
     782  if (Tess->LastTriangle != NULL) {
     783    // include the current position of the virtual sphere in the temporary raster3d file
     784    Vector *center = cloud->GetCenter();
     785    // make the circumsphere's center absolute again
     786    helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
     787    helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
     788    helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
     789    helper.Scale(1./3.);
     790    helper.SubtractVector(center);
     791    // and add to file plus translucency object
     792    *rasterfile << "# current virtual sphere\n";
     793    *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
     794    *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
     795    *rasterfile << "9\n  terminating special property\n";
     796    delete(center);
     797  }
    773798};
    774799
    … …  
    781806void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    782807{
     808        Info FunctionInfo(__func__);
    783809  TesselPoint *Walker = NULL;
    784810  int i;
    … …  
    826852void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N)
    827853{
     854        Info FunctionInfo(__func__);
    828855  if ((tecplot != NULL) && (TesselStruct != NULL)) {
    829856    // write header
    830857    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    831858    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
    832     *tecplot << "ZONE T=\"" << N << "-";
    833     for (int i=0;i<3;i++)
    834       *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
     859    *tecplot << "ZONE T=\"";
     860    if (N < 0) {
     861      *tecplot << cloud->GetName();
     862    } else {
     863      *tecplot << N << "-";
     864      for (int i=0;i<3;i++)
     865        *tecplot << (i==0 ? "" : "_") << TesselStruct->LastTriangle->endpoints[i]->node->Name;
     866    }
    835867    *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    836868    int i=0;
    … …  
    841873
    842874    // print atom coordinates
    843     Log() << Verbose(2) << "The following triangles were created:";
    844875    int Counter = 1;
    845876    TesselPoint *Walker = NULL;
    … …  
    851882    *tecplot << endl;
    852883    // print connectivity
     884    Log() << Verbose(1) << "The following triangles were created:" << endl;
    853885    for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
    854       Log() << Verbose(0) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name;
     886      Log() << Verbose(1) << " " << runner->second->endpoints[0]->node->Name << "<->" << runner->second->endpoints[1]->node->Name << "<->" << runner->second->endpoints[2]->node->Name << endl;
    855887      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    856888    }
    857889    delete[] (LookupList);
    858     Log() << Verbose(0) << endl;
    859890  }
    860891};
    … …  
    867898void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct)
    868899{
     900        Info FunctionInfo(__func__);
    869901  class BoundaryPointSet *point = NULL;
    870902  class BoundaryLineSet *line = NULL;
    871903
    872   //Log() << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
    873904  // calculate remaining concavity
    874905  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
    … …  
    878909    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    879910      line = LineRunner->second;
    880       //Log() << Verbose(2) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
     911      //Log() << Verbose(1) << "INFO: Current line of point " << *point << " is " << *line << "." << endl;
    881912      if (!line->CheckConvexityCriterion())
    882913        point->value += 1;
    883914    }
    884915  }
    885   //Log() << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;
    886916};
    887917
    … …  
    894924bool CheckListOfBaselines(const Tesselation * const TesselStruct)
    895925{
     926        Info FunctionInfo(__func__);
    896927  LineMap::const_iterator testline;
    897928  bool result = false;
    … …  
    901932  for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
    902933    if (testline->second->triangles.size() != 2) {
    903       Log() << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl;
     934      Log() << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl;
    904935      counter++;
    905936    }
    … …  
    912943}
    913944
     945/** Counts the number of triangle pairs that contain the given polygon.
     946 * \param *P polygon with endpoints to look for
     947 * \param *T set of triangles to create pairs from containing \a *P
     948 */
     949int CountTrianglePairContainingPolygon(const BoundaryPolygonSet * const P, const TriangleSet * const T)
     950{
     951  Info FunctionInfo(__func__);
     952  // check number of endpoints in *P
     953  if (P->endpoints.size() != 4) {
     954    eLog() << Verbose(1) << "CountTrianglePairContainingPolygon works only on polygons with 4 nodes!" << endl;
     955    return 0;
     956  }
     957
     958  // check number of triangles in *T
     959  if (T->size() < 2) {
     960    eLog() << Verbose(1) << "Not enough triangles to have pairs!" << endl;
     961    return 0;
     962  }
     963
     964  Log() << Verbose(0) << "Polygon is " << *P << endl;
     965  // create each pair, get the endpoints and check whether *P is contained.
     966  int counter = 0;
     967  PointSet Trianglenodes;
     968  class BoundaryPolygonSet PairTrianglenodes;
     969  for(TriangleSet::iterator Walker = T->begin(); Walker != T->end(); Walker++) {
     970    for (int i=0;i<3;i++)
     971      Trianglenodes.insert((*Walker)->endpoints[i]);
     972
     973    for(TriangleSet::iterator PairWalker = Walker; PairWalker != T->end(); PairWalker++) {
     974      if (Walker != PairWalker) { // skip first
     975        PairTrianglenodes.endpoints = Trianglenodes;
     976        for (int i=0;i<3;i++)
     977          PairTrianglenodes.endpoints.insert((*PairWalker)->endpoints[i]);
     978        const int size = PairTrianglenodes.endpoints.size();
     979        if (size == 4) {
     980          Log() << Verbose(0) << " Current pair of triangles: " << **Walker << "," << **PairWalker << " with " << size << " distinct endpoints:" << PairTrianglenodes << endl;
     981          // now check
     982          if (PairTrianglenodes.ContainsPresentTupel(P)) {
     983            counter++;
     984            Log() << Verbose(0) << "  ACCEPT: Matches with " << *P << endl;
     985          } else {
     986            Log() << Verbose(0) << "  REJECT: No match with " << *P << endl;
     987          }
     988        } else {
     989          Log() << Verbose(0) << "  REJECT: Less than four endpoints." << endl;
     990        }
     991      }
     992    }
     993    Trianglenodes.clear();
     994  }
     995  return counter;
     996};
     997
     998/** Checks whether two give polygons have two or more points in common.
     999 * \param *P1 first polygon
     1000 * \param *P2 second polygon
     1001 * \return true - are connected, false = are note
     1002 */
     1003bool ArePolygonsEdgeConnected(const BoundaryPolygonSet * const P1, const BoundaryPolygonSet * const P2)
     1004{
     1005  Info FunctionInfo(__func__);
     1006  int counter = 0;
     1007  for(PointSet::const_iterator Runner = P1->endpoints.begin(); Runner != P1->endpoints.end(); Runner++) {
     1008    if (P2->ContainsBoundaryPoint((*Runner))) {
     1009      counter++;
     1010      Log() << Verbose(1) << *(*Runner) << " of second polygon is found in the first one." << endl;
     1011      return true;
     1012    }
     1013  }
     1014  return false;
     1015};
     1016
     1017/** Combines second into the first and deletes the second.
     1018 * \param *P1 first polygon, contains all nodes on return
     1019 * \param *&P2 second polygon, is deleted.
     1020 */
     1021void CombinePolygons(BoundaryPolygonSet * const P1, BoundaryPolygonSet * &P2)
     1022{
     1023  Info FunctionInfo(__func__);
     1024  pair <PointSet::iterator, bool> Tester;
     1025  for(PointSet::iterator Runner = P2->endpoints.begin(); Runner != P2->endpoints.end(); Runner++) {
     1026    Tester = P1->endpoints.insert((*Runner));
     1027    if (Tester.second)
     1028      Log() << Verbose(0) << "Inserting endpoint " << *(*Runner) << " into first polygon." << endl;
     1029  }
     1030  P2->endpoints.clear();
     1031  delete(P2);
     1032};
     1033
Note: See TracChangeset for help on using the changeset viewer.