Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/tesselationhelpers.cpp

    rc5f836 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) {
    190194    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;
    … …  
    219223double GetPathLengthonCircumCircle(const Vector &CircleCenter, const Vector &CirclePlaneNormal, const double CircleRadius, const Vector &NewSphereCenter, const Vector &OldSphereCenter, const Vector &NormalVector, const Vector &SearchDirection)
    220224{
     225        Info FunctionInfo(__func__);
    221226  Vector helper;
    222227  double radius, alpha;
    223 
    224   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);
    225236  // test whether new center is on the parameter circle's plane
    226237  if (fabs(helper.ScalarProduct(&CirclePlaneNormal)) > HULLEPSILON) {
    … …  
    228239    helper.ProjectOntoPlane(&CirclePlaneNormal);
    229240  }
    230   radius = helper.ScalarProduct(&helper);
     241  radius = helper.NormSquared();
    231242  // test whether the new center vector has length of CircleRadius
    232243  if (fabs(radius - CircleRadius) > HULLEPSILON)
    233244    eLog() << Verbose(1) << "The projected center of the new sphere has radius " << radius << " instead of " << CircleRadius << "." << endl;
    234   alpha = helper.Angle(&OldSphereCenter);
     245  alpha = helper.Angle(&RelativeOldSphereCenter);
    235246  // make the angle unique by checking the halfplanes/search direction
    236247  if (helper.ScalarProduct(&SearchDirection) < -HULLEPSILON)  // acos is not unique on [0, 2.*M_PI), hence extra check to decide between two half intervals
    237248    alpha = 2.*M_PI - alpha;
    238   //Log() << Verbose(2) << "INFO: RelativeNewSphereCenter is " << helper << ", RelativeOldSphereCenter is " << OldSphereCenter << " and resulting angle is " << alpha << "." << endl;
    239   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);
    240251  helper.ProjectOntoPlane(&NormalVector);
    241252  // check whether new center is somewhat away or at least right over the current baseline to prevent intersecting triangles
    242253  if ((radius > HULLEPSILON) || (helper.Norm() < HULLEPSILON)) {
    243     //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;
    244255    return alpha;
    245256  } else {
    246     //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;
    247258    return 2.*M_PI;
    248259  }
    … …  
    264275double MinIntersectDistance(const gsl_vector * x, void *params)
    265276{
     277        Info FunctionInfo(__func__);
    266278  double retval = 0;
    267279  struct Intersection *I = (struct Intersection *)params;
    … …  
    284296
    285297  retval = HeightA.ScalarProduct(&HeightA) + HeightB.ScalarProduct(&HeightB);
    286   //Log() << Verbose(2) << "MinIntersectDistance called, result: " << retval << endl;
     298  //Log() << Verbose(1) << "MinIntersectDistance called, result: " << retval << endl;
    287299
    288300  return retval;
    … …  
    304316bool existsIntersection(const Vector &point1, const Vector &point2, const Vector &point3, const Vector &point4)
    305317{
     318        Info FunctionInfo(__func__);
    306319  bool result;
    307320
    … …  
    351364
    352365        if (status == GSL_SUCCESS) {
    353           Log() << Verbose(2) << "converged to minimum" <<  endl;
     366          Log() << Verbose(1) << "converged to minimum" <<  endl;
    354367        }
    355368    } while (status == GSL_CONTINUE && iter < 100);
    … …  
    376389  t2 = HeightB.ScalarProduct(&SideB)/SideB.ScalarProduct(&SideB);
    377390
    378   Log() << Verbose(2) << "Intersection " << intersection << " is at "
     391  Log() << Verbose(1) << "Intersection " << intersection << " is at "
    379392    << t1 << " for (" << point1 << "," << point2 << ") and at "
    380393    << t2 << " for (" << point3 << "," << point4 << "): ";
    381394
    382395  if (((t1 >= 0) && (t1 <= 1)) && ((t2 >= 0) && (t2 <= 1))) {
    383     Log() << Verbose(0) << "true intersection." << endl;
     396    Log() << Verbose(1) << "true intersection." << endl;
    384397    result = true;
    385398  } else {
    386     Log() << Verbose(0) << "intersection out of region of interest." << endl;
     399    Log() << Verbose(1) << "intersection out of region of interest." << endl;
    387400    result = false;
    388401  }
    … …  
    407420double GetAngle(const Vector &point, const Vector &reference, const Vector &OrthogonalVector)
    408421{
     422        Info FunctionInfo(__func__);
    409423  if (reference.IsZero())
    410424    return M_PI;
    … …  
    418432  }
    419433
    420   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;
    421435
    422436  return phi;
    … …  
    433447double CalculateVolumeofGeneralTetraeder(const Vector &a, const Vector &b, const Vector &c, const Vector &d)
    434448{
     449        Info FunctionInfo(__func__);
    435450  Vector Point, TetraederVector[3];
    436451  double volume;
    … …  
    456471bool CheckLineCriteriaForDegeneratedTriangle(const BoundaryPointSet * const nodes[3])
    457472{
     473        Info FunctionInfo(__func__);
    458474  bool result = false;
    459475  int counter = 0;
    … …  
    482498    }
    483499  if ((!result) && (counter > 1)) {
    484     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;
    485501    result = true;
    486502  }
    … …  
    489505
    490506
    491 /** Sort function for the candidate list.
    492  */
    493 bool SortCandidates(const CandidateForTesselation* candidate1, const CandidateForTesselation* candidate2)
    494 {
    495   Vector BaseLineVector, OrthogonalVector, helper;
    496   if (candidate1->BaseLine != candidate2->BaseLine) {  // sanity check
    497     eLog() << Verbose(1) << "sortCandidates was called for two different baselines: " << candidate1->BaseLine << " and " << candidate2->BaseLine << "." << endl;
    498     //return false;
    499     exit(1);
    500   }
    501   // create baseline vector
    502   BaseLineVector.CopyVector(candidate1->BaseLine->endpoints[1]->node->node);
    503   BaseLineVector.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    504   BaseLineVector.Normalize();
    505 
    506   // 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!)
    507   helper.CopyVector(candidate1->BaseLine->endpoints[0]->node->node);
    508   helper.SubtractVector(candidate1->point->node);
    509   OrthogonalVector.CopyVector(&helper);
    510   helper.VectorProduct(&BaseLineVector);
    511   OrthogonalVector.SubtractVector(&helper);
    512   OrthogonalVector.Normalize();
    513 
    514   // calculate both angles and correct with in-plane vector
    515   helper.CopyVector(candidate1->point->node);
    516   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    517   double phi = BaseLineVector.Angle(&helper);
    518   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    519     phi = 2.*M_PI - phi;
    520   }
    521   helper.CopyVector(candidate2->point->node);
    522   helper.SubtractVector(candidate1->BaseLine->endpoints[0]->node->node);
    523   double psi = BaseLineVector.Angle(&helper);
    524   if (OrthogonalVector.ScalarProduct(&helper) > 0) {
    525     psi = 2.*M_PI - psi;
    526   }
    527 
    528   Log() << Verbose(2) << *candidate1->point << " has angle " << phi << endl;
    529   Log() << Verbose(2) << *candidate2->point << " has angle " << psi << endl;
    530 
    531   // return comparison
    532   return phi < psi;
    533 };
     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//};
    534551
    535552/**
    … …  
    543560TesselPoint* FindSecondClosestPoint(const Vector* Point, const LinkedCell* const LC)
    544561{
     562        Info FunctionInfo(__func__);
    545563  TesselPoint* closestPoint = NULL;
    546564  TesselPoint* secondClosestPoint = NULL;
    … …  
    553571  for(int i=0;i<NDIM;i++) // store indices of this cell
    554572    N[i] = LC->n[i];
    555   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;
    556574
    557575  LC->GetNeighbourBounds(Nlower, Nupper);
    558   //Log() << Verbose(0) << endl;
     576  //Log() << Verbose(1) << endl;
    559577  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    560578    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    561579      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    562580        const LinkedNodes *List = LC->GetCurrentCell();
    563         //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;
    564582        if (List != NULL) {
    565583          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    … …  
    597615TesselPoint* FindClosestPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const LC)
    598616{
     617        Info FunctionInfo(__func__);
    599618  TesselPoint* closestPoint = NULL;
    600619  SecondPoint = NULL;
    … …  
    607626  for(int i=0;i<NDIM;i++) // store indices of this cell
    608627    N[i] = LC->n[i];
    609   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;
    610629
    611630  LC->GetNeighbourBounds(Nlower, Nupper);
    612   //Log() << Verbose(0) << endl;
     631  //Log() << Verbose(1) << endl;
    613632  for (LC->n[0] = Nlower[0]; LC->n[0] <= Nupper[0]; LC->n[0]++)
    614633    for (LC->n[1] = Nlower[1]; LC->n[1] <= Nupper[1]; LC->n[1]++)
    615634      for (LC->n[2] = Nlower[2]; LC->n[2] <= Nupper[2]; LC->n[2]++) {
    616635        const LinkedNodes *List = LC->GetCurrentCell();
    617         //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;
    618637        if (List != NULL) {
    619638          for (LinkedNodes::const_iterator Runner = List->begin(); Runner != List->end(); Runner++) {
    … …  
    626645              distance = currentNorm;
    627646              closestPoint = (*Runner);
    628               //Log() << Verbose(2) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
     647              //Log() << Verbose(1) << "INFO: New Nearest Neighbour is " << *closestPoint << "." << endl;
    629648            } else if (currentNorm < secondDistance) {
    630649              secondDistance = currentNorm;
    631650              SecondPoint = (*Runner);
    632               //Log() << Verbose(2) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
     651              //Log() << Verbose(1) << "INFO: New Second Nearest Neighbour is " << *SecondPoint << "." << endl;
    633652            }
    634653          }
    … …  
    640659  // output
    641660  if (closestPoint != NULL) {
    642     Log() << Verbose(2) << "Closest point is " << *closestPoint;
     661    Log() << Verbose(1) << "Closest point is " << *closestPoint;
    643662    if (SecondPoint != NULL)
    644663      Log() << Verbose(0) << " and second closest is " << *SecondPoint;
    … …  
    656675Vector * GetClosestPointBetweenLine(const BoundaryLineSet * const Base, const BoundaryLineSet * const OtherBase)
    657676{
     677        Info FunctionInfo(__func__);
    658678  // construct the plane of the two baselines (i.e. take both their directional vectors)
    659679  Vector Normal;
    … …  
    666686  Normal.VectorProduct(&OtherBaseline);
    667687  Normal.Normalize();
    668   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;
    669689
    670690  // project one offset point of OtherBase onto this plane (and add plane offset vector)
    … …  
    683703  Normal.CopyVector(Intersection);
    684704  Normal.SubtractVector(Base->endpoints[0]->node->node);
    685   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;
    686706
    687707  return Intersection;
    … …  
    696716double DistanceToTrianglePlane(const Vector *x, const BoundaryTriangleSet * const triangle)
    697717{
     718        Info FunctionInfo(__func__);
    698719  double distance = 0.;
    699720  if (x == NULL) {
    … …  
    712733void WriteVrmlFile(ofstream * const vrmlfile, const Tesselation * const Tess, const PointCloud * const cloud)
    713734{
     735        Info FunctionInfo(__func__);
    714736  TesselPoint *Walker = NULL;
    715737  int i;
    … …  
    755777void IncludeSphereinRaster3D(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    756778{
     779        Info FunctionInfo(__func__);
    757780  Vector helper;
    758   // include the current position of the virtual sphere in the temporary raster3d file
    759   Vector *center = cloud->GetCenter();
    760   // make the circumsphere's center absolute again
    761   helper.CopyVector(Tess->LastTriangle->endpoints[0]->node->node);
    762   helper.AddVector(Tess->LastTriangle->endpoints[1]->node->node);
    763   helper.AddVector(Tess->LastTriangle->endpoints[2]->node->node);
    764   helper.Scale(1./3.);
    765   helper.SubtractVector(center);
    766   // and add to file plus translucency object
    767   *rasterfile << "# current virtual sphere\n";
    768   *rasterfile << "8\n  25.0    0.6     -1.0 -1.0 -1.0     0.2        0 0 0 0\n";
    769   *rasterfile << "2\n  " << helper.x[0] << " " << helper.x[1] << " " << helper.x[2] << "\t" << 5. << "\t1 0 0\n";
    770   *rasterfile << "9\n  terminating special property\n";
    771   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  }
    772798};
    773799
    … …  
    780806void WriteRaster3dFile(ofstream * const rasterfile, const Tesselation * const Tess, const PointCloud * const cloud)
    781807{
     808        Info FunctionInfo(__func__);
    782809  TesselPoint *Walker = NULL;
    783810  int i;
    … …  
    825852void WriteTecplotFile(ofstream * const tecplot, const Tesselation * const TesselStruct, const PointCloud * const cloud, const int N)
    826853{
     854        Info FunctionInfo(__func__);
    827855  if ((tecplot != NULL) && (TesselStruct != NULL)) {
    828856    // write header
    829857    *tecplot << "TITLE = \"3D CONVEX SHELL\"" << endl;
    830858    *tecplot << "VARIABLES = \"X\" \"Y\" \"Z\" \"U\"" << endl;
    831     *tecplot << "ZONE T=\"" << N << "-";
    832     for (int i=0;i<3;i++)
    833       *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    }
    834867    *tecplot << "\", N=" << TesselStruct->PointsOnBoundary.size() << ", E=" << TesselStruct->TrianglesOnBoundary.size() << ", DATAPACKING=POINT, ZONETYPE=FETRIANGLE" << endl;
    835868    int i=0;
    … …  
    840873
    841874    // print atom coordinates
    842     Log() << Verbose(2) << "The following triangles were created:";
    843875    int Counter = 1;
    844876    TesselPoint *Walker = NULL;
    … …  
    850882    *tecplot << endl;
    851883    // print connectivity
     884    Log() << Verbose(1) << "The following triangles were created:" << endl;
    852885    for (TriangleMap::const_iterator runner = TesselStruct->TrianglesOnBoundary.begin(); runner != TesselStruct->TrianglesOnBoundary.end(); runner++) {
    853       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;
    854887      *tecplot << LookupList[runner->second->endpoints[0]->node->nr] << " " << LookupList[runner->second->endpoints[1]->node->nr] << " " << LookupList[runner->second->endpoints[2]->node->nr] << endl;
    855888    }
    856889    delete[] (LookupList);
    857     Log() << Verbose(0) << endl;
    858890  }
    859891};
    … …  
    866898void CalculateConcavityPerBoundaryPoint(const Tesselation * const TesselStruct)
    867899{
     900        Info FunctionInfo(__func__);
    868901  class BoundaryPointSet *point = NULL;
    869902  class BoundaryLineSet *line = NULL;
    870903
    871   //Log() << Verbose(2) << "Begin of CalculateConcavityPerBoundaryPoint" << endl;
    872904  // calculate remaining concavity
    873905  for (PointMap::const_iterator PointRunner = TesselStruct->PointsOnBoundary.begin(); PointRunner != TesselStruct->PointsOnBoundary.end(); PointRunner++) {
    … …  
    877909    for (LineMap::iterator LineRunner = point->lines.begin(); LineRunner != point->lines.end(); LineRunner++) {
    878910      line = LineRunner->second;
    879       //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;
    880912      if (!line->CheckConvexityCriterion())
    881913        point->value += 1;
    882914    }
    883915  }
    884   //Log() << Verbose(2) << "End of CalculateConcavityPerBoundaryPoint" << endl;
    885916};
    886917
    … …  
    893924bool CheckListOfBaselines(const Tesselation * const TesselStruct)
    894925{
     926        Info FunctionInfo(__func__);
    895927  LineMap::const_iterator testline;
    896928  bool result = false;
    … …  
    900932  for (testline = TesselStruct->LinesOnBoundary.begin(); testline != TesselStruct->LinesOnBoundary.end(); testline++) {
    901933    if (testline->second->triangles.size() != 2) {
    902       Log() << Verbose(1) << *testline->second << "\t" << testline->second->triangles.size() << endl;
     934      Log() << Verbose(2) << *testline->second << "\t" << testline->second->triangles.size() << endl;
    903935      counter++;
    904936    }
    … …  
    911943}
    912944
     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.