source: molecuilder/src/tesselation.hpp@ ff4611

/*
 * 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;

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

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

#define ParallelEpsilon 1e-3

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

#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 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 *>

/********************************************** 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 * Walker);
    ~BoundaryPointSet();

    void AddLine(class BoundaryLineSet *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(class BoundaryPointSet *Point[2], const int number);
    ~BoundaryLineSet();

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

    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 *line[3], int number);
    ~BoundaryTriangleSet();

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

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

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 TesselPoint *GetTerminalPoint() const { return NULL; };
  virtual int GetMaxId() const { return 0; };
  virtual void GoToNext() const {};
  virtual void GoToPrevious() const {};
  virtual void GoToFirst() const {};
  virtual void GoToLast() 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, Vector OptCandidateCenter, Vector OtherOptCandidateCenter);
  ~CandidateForTesselation();

  TesselPointList pointlist;
  BoundaryLineSet *BaseLine;
  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(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
    void AlwaysAddTesselationTriangleLine(class BoundaryPointSet *a, class BoundaryPointSet *b, const int n);
    void AddTesselationTriangle();
    void AddTesselationTriangle(const int nr);
    void AddCandidateTriangle(CandidateForTesselation CandidateLine);
    void RemoveTesselationTriangle(class BoundaryTriangleSet *triangle);
    void RemoveTesselationLine(class BoundaryLineSet *line);
    void RemoveTesselationPoint(class BoundaryPointSet *point);


    // concave envelope
    void 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(Vector &NormalVector, Vector &SearchDirection, Vector &OldSphereCenter, CandidateForTesselation &CandidateLine, const class TesselPoint  * const ThirdNode, const double RADIUS, const LinkedCell *LC) const;
    bool FindNextSuitableTriangle(CandidateForTesselation &CandidateLine, BoundaryTriangleSet &T, const double& RADIUS, const LinkedCell *LC);
    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);
    map<int, int> * FindAllDegeneratedTriangles();
    map<int, int> * FindAllDegeneratedLines();
    void RemoveDegeneratedTriangles();
    void AddBoundaryPointByDegeneratedTriangle(class TesselPoint *point, LinkedCell *LC);
    int CorrectAllDegeneratedPolygons();

    set<TesselPoint*> * GetAllConnectedPoints(const TesselPoint* const Point) const;
    set<BoundaryTriangleSet*> *GetAllTriangles(const BoundaryPointSet * const Point) const;
    list<list<TesselPoint*> *> * GetPathsOfConnectedPoints(const TesselPoint* const Point) const;
    list<list<TesselPoint*> *> * GetClosedPathsOfConnectedPoints(const TesselPoint* const Point) const;
    list<TesselPoint*> * GetCircleOfSetOfPoints(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference = NULL) const;
    list<TesselPoint*> * GetCircleOfConnectedTriangles(set<TesselPoint*> *SetOfNeighbours, const TesselPoint* const Point, const Vector * const Reference) const;
    class BoundaryPointSet *GetCommonEndpoint(const BoundaryLineSet * line1, const BoundaryLineSet * line2) const;
    list<BoundaryTriangleSet*> *FindTriangles(const TesselPoint* const Points[3]) const;
    list<BoundaryTriangleSet*> * FindClosestTrianglesToPoint(const Vector *x, const LinkedCell* LC) const;
    class BoundaryTriangleSet * FindClosestTriangleToPoint(const Vector *x, const LinkedCell* LC) const;
    bool IsInnerPoint(const Vector &Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
    bool IsInnerPoint(const TesselPoint * const Point, const LinkedCell* const LC, const double epsilon = -MYEPSILON) const;
    bool AddBoundaryPoint(TesselPoint * Walker, const int n);
    TesselPoint* FindSecondClosestBoundaryPoint(const Vector* Point, const LinkedCell* const LC) const;
    TesselPoint* FindClosestBoundaryPoint(const Vector* Point, TesselPoint *&SecondPoint, const LinkedCell* const 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 TesselPoint *GetTerminalPoint() const;
    virtual void GoToNext() const;
    virtual void GoToPrevious() const;
    virtual void GoToFirst() const;
    virtual void GoToLast() 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;
};


#endif /* TESSELATION_HPP_ */