source: src/tesselation.hpp@ 88104f

/*
 * tesselation.hpp
 *
 *  The tesselation class is meant to contain the envelope (concave, convex or neither) of a set of Vectors. As we actually mean this stuff for atoms, we have to encapsulate it all a bit.
 *
 *  Created on: Aug 3, 2009
 *      Author: heber
 */

#ifndef TESSELATION_HPP_
#define TESSELATION_HPP_

using namespace std;

/*********************************************** includes ***********************************/

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <map>
#include <list>
#include <set>
#include <stack>

#include "atom_particleinfo.hpp"
#include "helpers.hpp"
#include "vector.hpp"

/****************************************** forward declarations *****************************/

class BoundaryPointSet;
class BoundaryLineSet;
class BoundaryTriangleSet;
class LinkedCell;
class TesselPoint;
class PointCloud;
class Tesselation;
class Plane;

/********************************************** definitions *********************************/

#define DoTecplotOutput 1
#define DoRaster3DOutput 1
#define DoVRMLOutput 0
#define TecplotSuffix ".dat"
#define Raster3DSuffix ".r3d"
#define VRMLSUffix ".wrl"

#define ParallelEpsilon 1e-3

// ======================================================= some template functions =========================================

#define IndexToIndex map <int, int>

#define PointMap map < int, class BoundaryPointSet * >
#define PointSet set < class BoundaryPointSet * >
#define PointList list < class BoundaryPointSet * >
#define PointPair pair < int, class BoundaryPointSet * >
#define PointTestPair pair < PointMap::iterator, bool >

#define CandidateList list <class CandidateForTesselation *>
#define CandidateMap map <class BoundaryLineSet *, class CandidateForTesselation *>

#define LineMap multimap < int, class BoundaryLineSet * >
#define LineSet set < class BoundaryLineSet * >
#define LineList list < class BoundaryLineSet * >
#define LinePair pair < int, class BoundaryLineSet * >
#define LineTestPair pair < LineMap::iterator, bool >

#define TriangleMap map < int, class BoundaryTriangleSet * >
#define TriangleSet set < class BoundaryTriangleSet * >
#define TriangleList list < class BoundaryTriangleSet * >
#define TrianglePair pair < int, class BoundaryTriangleSet * >
#define TriangleTestPair pair < TrianglePair::iterator, bool >

#define PolygonMap map < int, class BoundaryPolygonSet * >
#define PolygonSet set < class BoundaryPolygonSet * >
#define PolygonList list < class BoundaryPolygonSet * >

#define DistanceToPointMap multimap <double, class BoundaryPointSet * >
#define DistanceToPointPair pair <double, class BoundaryPointSet * >

#define DistanceMultiMap multimap <double, pair < PointMap::iterator, PointMap::iterator> >
#define DistanceMultiMapPair pair <double, pair < PointMap::iterator, PointMap::iterator> >

#define TesselPointList list <TesselPoint *>
#define TesselPointSet set <TesselPoint *>

#define ListOfTesselPointList list<list <TesselPoint *> *>

enum centers {Opt, OtherOpt};

/********************************************** declarations *******************************/

template <typename T> void SetEndpointsOrdered(T endpoints[2], T endpoint1, T endpoint2)
{
  if (endpoint1->Nr < endpoint2->Nr) {
    endpoints[0] = endpoint1;
    endpoints[1] = endpoint2;
  } else {
    endpoints[0] = endpoint2;
    endpoints[1] = endpoint1;
  }
};

// ======================================================== class BoundaryPointSet =========================================

class BoundaryPointSet {
  public:
    BoundaryPointSet();
    BoundaryPointSet(TesselPoint * const Walker);
    ~BoundaryPointSet();

    void AddLine(BoundaryLineSet * const line);

    LineMap lines;
    int LinesCount;
    TesselPoint *node;
    double value;
    int Nr;
};

ostream & operator << (ostream &ost, const BoundaryPointSet &a);

// ======================================================== class BoundaryLineSet ==========================================

class BoundaryLineSet {
  public:
    BoundaryLineSet();
    BoundaryLineSet(BoundaryPointSet * const Point[2], const int number);
    BoundaryLineSet(BoundaryPointSet * const Point1, BoundaryPointSet * const Point2, const int number);
    ~BoundaryLineSet();

    void AddTriangle(BoundaryTriangleSet * const triangle);
    bool IsConnectedTo(const BoundaryLineSet * const line) const;
    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
    bool CheckConvexityCriterion() const;
    class BoundaryPointSet *GetOtherEndpoint(const BoundaryPointSet * const point) const;

    class BoundaryPointSet *endpoints[2];
    TriangleMap triangles;
    int Nr;
    bool skipped;
};

ostream & operator << (ostream &ost, const BoundaryLineSet &a);

// ======================================================== class BoundaryTriangleSet =======================================

class BoundaryTriangleSet {
  public:
    BoundaryTriangleSet();
    BoundaryTriangleSet(class BoundaryLineSet * const line[3], const int number);
    ~BoundaryTriangleSet();

    void GetNormalVector(const Vector &NormalVector);
    void GetCenter(Vector * const center) const;
    bool GetIntersectionInsideTriangle(const Vector * const MolCenter, const Vector * const x, Vector * const Intersection) const;
    double GetClosestPointInsideTriangle(const Vector * const x, Vector * const ClosestPoint) const;
    bool ContainsBoundaryLine(const BoundaryLineSet * const line) const;
    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
    bool ContainsBoundaryPoint(const TesselPoint * const point) const;
    class BoundaryPointSet *GetThirdEndpoint(const BoundaryLineSet * const line) const;
    bool IsPresentTupel(const BoundaryPointSet * const Points[3]) const;
    bool IsPresentTupel(const BoundaryTriangleSet * const T) const;

    Plane getPlane() const;
    Vector getEndpoint(int) const;
    std::string getEndpointName(int) const;

    class BoundaryPointSet *endpoints[3];
    class BoundaryLineSet *lines[3];
    Vector NormalVector;
    Vector SphereCenter;
    int Nr;

  private:

};

ostream & operator << (ostream &ost, const BoundaryTriangleSet &a);


// ======================================================== class BoundaryTriangleSet =======================================

/** Set of BoundaryPointSet.
 * This is just meant as a container for a group of endpoints, extending the node, line, triangle concept. However, this has
 * only marginally something to do with the tesselation. Hence, there is no incorporation into the bookkeeping of the Tesselation
 * class (i.e. no allocation, no deletion).
 * \note we assume that the set of endpoints reside (more or less) on a plane.
 */
class BoundaryPolygonSet {
  public:
    BoundaryPolygonSet();
    ~BoundaryPolygonSet();

    Vector * GetNormalVector(const Vector &NormalVector) const;
    void GetCenter(Vector *center) const;
    bool ContainsBoundaryLine(const BoundaryLineSet * const line) const;
    bool ContainsBoundaryPoint(const BoundaryPointSet * const point) const;
    bool ContainsBoundaryPoint(const TesselPoint * const point) const;
    bool ContainsBoundaryTriangle(const BoundaryTriangleSet * const point) const;
    bool ContainsPresentTupel(const BoundaryPointSet * const * Points, const int dim) const;
    bool ContainsPresentTupel(const BoundaryPolygonSet * const P) const;
    bool ContainsPresentTupel(const PointSet &endpoints) const;
    TriangleSet * GetAllContainedTrianglesFromEndpoints() const;
    bool FillPolygonFromTrianglesOfLine(const BoundaryLineSet * const line);

    PointSet endpoints;
    int Nr;
};

ostream & operator << (ostream &ost, const BoundaryPolygonSet &a);

// =========================================================== class TESSELPOINT ===========================================

/** Is a single point of the set of Vectors, also a super-class to be inherited and and its functions to be implemented.
 */
class TesselPoint : virtual public ParticleInfo {
public:
  TesselPoint();
  virtual ~TesselPoint();

  Vector *node;   // pointer to position of the dot in space

  virtual ostream & operator << (ostream &ost);
};

ostream & operator << (ostream &ost, const TesselPoint &a);

// =========================================================== class POINTCLOUD ============================================

/** Super-class for all point clouds structures, also molecules. They have to inherit this structure and implement the virtual function to access the Vectors.
 * This basically encapsulates a list structure.
 */
class PointCloud {
public:
  PointCloud();
  virtual ~PointCloud();

  virtual const char * const GetName() const { return "unknown"; };
  virtual Vector *GetCenter() const { return NULL; };
  virtual TesselPoint *GetPoint() const { return NULL; };
  virtual int GetMaxId() const { return 0; };
  virtual void GoToNext() const {};
  virtual void GoToFirst() const {};
  virtual bool IsEmpty() const { return true; };
  virtual bool IsEnd() const { return true; };
};

// ======================================================== class CandidateForTesselation =========================================

class CandidateForTesselation {
  public :
  CandidateForTesselation(BoundaryLineSet* currentBaseLine);
  CandidateForTesselation(TesselPoint* candidate, BoundaryLineSet* currentBaseLine, BoundaryPointSet *point, Vector OptCandidateCenter, Vector OtherOptCandidateCenter);
  ~CandidateForTesselation();

  bool CheckValidity(const double RADIUS, const LinkedCell *LC) const;

  TesselPointList pointlist;
  const BoundaryLineSet * BaseLine;
  const BoundaryPointSet * ThirdPoint;
  const BoundaryTriangleSet *T;
  Vector OldCenter;
  Vector OptCenter;
  Vector OtherOptCenter;
  double ShortestAngle;
  double OtherShortestAngle;
};

ostream & operator <<(ostream &ost, const  CandidateForTesselation &a);

// =========================================================== class TESSELATION ===========================================

/** Contains the envelope to a PointCloud.
 */
class Tesselation : public PointCloud {
  public:

    Tesselation();
    virtual ~Tesselation();

    void AddTesselationPoint(TesselPoint* Candidate, const int n);
    void SetTesselationPoint(TesselPoint* Candidate, const int n) const;
    void AddTesselationLine(const Vector * const OptCenter, const BoundaryPointSet * const candidate, class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
    void AddNewTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
    void AddExistingTesselationTriangleLine(class BoundaryLineSet *FindLine, int n);
    void AddTesselationTriangle();
    void AddTesselationTriangle(const int nr);
    void AddCandidateTriangle(CandidateForTesselation &CandidateLine, enum centers type);
    void AddDegeneratedTriangle(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC);
    void AddCandidatePolygon(CandidateForTesselation CandidateLine, const double RADIUS, const LinkedCell *LC);
    void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle);
    void RemoveTesselationLine(class BoundaryLineSet *line);
    void RemoveTesselationPoint(class BoundaryPointSet *point);
    bool CheckDegeneracy(CandidateForTesselation &CandidateLine, const double RADIUS, const LinkedCell *LC) const;


    // concave envelope
    bool FindStartingTriangle(const double RADIUS, const LinkedCell *LC);
    void FindSecondPointForTesselation(class TesselPoint* a, Vector Oben, class TesselPoint*& OptCandidate, double Storage[3], double RADIUS, const LinkedCell *LC);
    void FindThirdPointForTesselation(const Vector &NormalVector, const Vector &SearchDirection, const Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class BoundaryPointSet  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const;
    bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, const BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
    bool FindCandidatesforOpenLines(const double RADIUS, const LinkedCell *&LCList);
    int CheckPresenceOfTriangle(class TesselPoint *Candidates[3]) const;
    class BoundaryTriangleSet * GetPresentTriangle(TesselPoint *Candidates[3]);

    // convex envelope
    void TesselateOnBoundary(const PointCloud * const cloud);
    void GuessStartingTriangle();
    bool InsertStraddlingPoints(const PointCloud *cloud, const LinkedCell *LC);
    double RemovePointFromTesselatedSurface(class BoundaryPointSet *point);
    class BoundaryLineSet * FlipBaseline(class BoundaryLineSet *Base);
    double PickFarthestofTwoBaselines(class BoundaryLineSet *Base);
    class BoundaryPointSet *IsConvexRectangle(class BoundaryLineSet *Base);
    IndexToIndex * FindAllDegeneratedTriangles();
    IndexToIndex * FindAllDegeneratedLines();
    void RemoveDegeneratedTriangles();
    void AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC);
    int CorrectAllDegeneratedPolygons();

    TesselPointSet * GetAllConnectedPoints(const TesselPoint* const Point) const;
    TriangleSet * GetAllTriangles(const BoundaryPointSet * const Point) const;
    ListOfTesselPointList * GetPathsOfConnectedPoints(const TesselPoint* const Point) const;
    ListOfTesselPointList * GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const;
    TesselPointList * GetCircleOfSetOfPoints(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
    TesselPointList * GetCircleOfConnectedTriangles(TesselPointSet *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const;
    class BoundaryPointSet * GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
    TriangleList * FindTriangles(const TesselPoint* const Points[3]) const;
    TriangleList * FindClosestTrianglesToVector(const Vector *x, const LinkedCell* LC) const;
    BoundaryTriangleSet * FindClosestTriangleToVector(const Vector *x, const LinkedCell* LC) const;
    bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC) const;
    double GetDistanceSquaredToTriangle(const Vector &Point, const BoundaryTriangleSet* const triangle) const;
    double GetDistanceToSurface(const Vector &Point, const LinkedCell* const LC) const;
    BoundaryTriangleSet * GetClosestTriangleOnSurface(const Vector &Point, const LinkedCell* const LC) const;
    bool AddBoundaryPoint(TesselPoint * Walker, const int n);
    DistanceToPointMap * FindClosestBoundaryPointsToVector(const Vector *x, const LinkedCell* LC) const;
    BoundaryLineSet * FindClosestBoundaryLineToVector(const Vector *x, const LinkedCell* LC) const;

    // print for debugging
    void PrintAllBoundaryPoints(ofstream *out) const;
    void PrintAllBoundaryLines(ofstream *out) const;
    void PrintAllBoundaryTriangles(ofstream *out) const;

    // store envelope in file
    void Output(const char *filename, const PointCloud * const cloud);

    PointMap PointsOnBoundary;
    LineMap LinesOnBoundary;
    CandidateMap OpenLines;
    TriangleMap TrianglesOnBoundary;
    int PointsOnBoundaryCount;
    int LinesOnBoundaryCount;
    int TrianglesOnBoundaryCount;

    // PointCloud implementation for PointsOnBoundary
    virtual Vector *GetCenter(ofstream *out) const;
    virtual TesselPoint *GetPoint() const;
    virtual void GoToNext() const;
    virtual void GoToFirst() const;
    virtual bool IsEmpty() const;
    virtual bool IsEnd() const;

    class BoundaryPointSet *BPS[2];
    class BoundaryLineSet *BLS[3];
    class BoundaryTriangleSet *BTS;
    class BoundaryTriangleSet *LastTriangle;
    int TriangleFilesWritten;

  private:
    mutable class BoundaryPointSet *TPS[3]; //this is a Storage for pointers to triangle points, this and BPS[2] needed due to AddLine restrictions

    mutable PointMap::const_iterator InternalPointer;

    //bool HasOtherBaselineBetterCandidate(const BoundaryLineSet * const BaseRay, const TesselPoint * const OptCandidate, double ShortestAngle, double RADIUS, const LinkedCell * const LC) const;
    void FindDegeneratedCandidatesforOpenLines(TesselPoint * const Sprinter, const Vector * const OptCenter);
};


#endif /* TESSELATION_HPP_ */