delaunay.h

Go to the documentation of this file.

#ifndef delaunay_h
#define delaunay_h

/*+
________________________________________________________________________

 (C) dGB Beheer B.V.; (LICENSE) http://opendtect.org/OpendTect_license.txt
 Author:        Y.C. Liu
 Date:          January 2008
 RCS:           $Id$
________________________________________________________________________

-*/

#include "algomod.h"
#include "coord.h"
#include "odmemory.h"
#include "typeset.h"
#include "thread.h"
#include "trigonometry.h"
class od_ostream;


#define mDAGTriangleForceSingleThread

mExpClass(Algo) DAGTriangleTree
{
public:
                        DAGTriangleTree();
                        DAGTriangleTree(const DAGTriangleTree&);
    virtual             ~DAGTriangleTree();
    DAGTriangleTree&    operator=(const DAGTriangleTree&);

    static bool         computeCoordRanges(const TypeSet<Coord>&,
                            Interval<double>&,Interval<double>&);

    bool                setCoordList(const TypeSet<Coord>* OD);
    bool                setCoordList(TypeSet<Coord>*,OD::PtrPolicy);
    const TypeSet<Coord>& coordList() const { return *coordlist_; }

    bool                setBBox(const Interval<double>& xrg,
                                const Interval<double>& yrg);

    bool                isOK() const { return triangles_.size(); }

    bool                init();

    bool                insertPoint(int pointidx, int& dupid);
    int                 insertPoint(const Coord&, int& dupid);

    const Coord         getInitCoord(int vetexidx) const;
    bool                getTriangle(const Coord&,int& dupid,
                                    TypeSet<int>& vertexindices) const;
    bool                getCoordIndices(TypeSet<int>&) const;
    bool                getSurroundingIndices(TypeSet<int>&) const;
    bool                getConnections(int pointidx,TypeSet<int>&) const;
    bool                getWeights(int pointidx,const TypeSet<int>& conns,
                                   TypeSet<double>& weights,
                                   bool normalize=true) const;
    bool                getConnectionAndWeights(int ptidx,TypeSet<int>& conns,
                                             TypeSet<double>& weights,
                                             bool normailze=true) const;
    void                setEpsilon(double err)  { epsilon_ = err; }

    void                dumpTo(od_ostream&) const;
    void                dumpTriangulationToIV(od_ostream&) const;


    static int          cNoVertex()     { return -1; }

protected:

    static char         cIsOutside()    { return 0; }
    static char         cIsInside()     { return 1; }
    static char         cIsDuplicate()  { return 2; }
    static char         cError()        { return -1; }

    static int          cNoTriangle()   { return -1; }
    static int          cInitVertex0()  { return -2; }
    static int          cInitVertex1()  { return -3; }
    static int          cInitVertex2()  { return -4; }

    char        searchTriangle(const Coord& pt,int start,int& t0,
                               int& dupid) const;
    char        searchFurther(const Coord& pt,int& nti0,int& dupid) const;

    void        splitTriangleInside(int ci,int ti);
    void        legalizeTriangles(TypeSet<char>& v0s,TypeSet<char>& v1s,
                        TypeSet<int>& tis);
    int         getNeighbor(int v0,int v1,int ti) const;
    int         searchChild(int v0,int v1,int ti) const;
    char        isInside(const Coord& pt,int ti,int& dupid) const;

    struct DAGTriangle
    {
                        DAGTriangle();
        bool            operator==(const DAGTriangle&) const;
        DAGTriangle&    operator=(const DAGTriangle&);

        int             coordindices_[3];
        int             childindices_[3];
        int             neighbors_[3];

        bool            hasChildren() const;
    };

    bool                                multithreadsupport_;

    mutable Threads::ReadWriteLock      trianglelock_;

    double                              epsilon_;
    TypeSet<DAGTriangle>                triangles_;

    TypeSet<Coord>*                     coordlist_;
    bool                                ownscoordlist_;

    mutable Threads::ReadWriteLock      coordlock_;

    Coord                               initialcoords_[3];
};


mExpClass(Algo) DelaunayTriangulator : public ParallelTask
{ mODTextTranslationClass(DelaunayTriangulator);
public:
                        DelaunayTriangulator(DAGTriangleTree&);
                        ~DelaunayTriangulator();

    bool                isDataRandom()          { return israndom_; }
    void                dataIsRandom(bool yn)   { israndom_ = yn; }
    void                setCalcScope(const Interval<int>& rg);

protected:

#ifdef mDAGTriangleForceSingleThread
    int                 maxNrThreads() const { return 1; }
#endif
    od_int64            nrIterations() const;
    bool                doWork(od_int64,od_int64, int );
    bool                doPrepare(int);

    uiString            uiNrDoneText() const
                        { return tr("Points triangulated"); }
    uiString            uiMessage() const { return tr("Triangulating"); }

    od_int64*           permutation_;
    bool                israndom_;
    DAGTriangleTree&    tree_;
    Interval<od_int64>  calcscope_;
};


mExpClass(Algo) Triangle2DInterpolator
{
public:
                        Triangle2DInterpolator(const DAGTriangleTree&);

                        /*The vertices are indices from the DAGTriangleTree
                          coordlist, corresponding to the weights.*/
    bool                computeWeights(const Coord&,TypeSet<int>& vertices,
                                       TypeSet<double>& weights,
                                       double maxdist=mUdf(double),
                                       bool dointerpolate=true);
                        /*If don't do interpolate, nearest node will be taken.*/

protected:

    bool                        setFromAzimuth(const TypeSet<int>& tmpvertices,
                                               const Coord&,
                                               TypeSet<int>& vertices,
                                               TypeSet<double>& weights);
    const DAGTriangleTree&      triangles_;
    TypeSet<int>                corner0_;
    TypeSet<double>             cornerweights0_;
    TypeSet<int>                corner1_;
    TypeSet<double>             cornerweights1_;
    TypeSet<int>                corner2_;
    TypeSet<double>             cornerweights2_;
    Coord                       initcenter_;

    TypeSet<int>                perimeter_;
    TypeSet<double>             perimeterazimuth_;
};


/*Simple polyon triangulation, does not work if you have holes inside it.
  return each triangle with three indicies in order. */
inline  bool PolygonTriangulate( const TypeSet<Coord>& knots,TypeSet<int>& res )
{
    const int nrknots = knots.size();
    if ( nrknots < 3 )
      return false;

    /* Make sure it is a counter-clockwise polygon in ci */
    double area=0;
    for( int idx=nrknots-1, idy=0; idy<nrknots; idx=idy++ )
      area += (knots[idx].x*knots[idy].y - knots[idy].x*knots[idx].y);
    area *= 0.5;

    TypeSet<int> ci;
    if ( 0<area )
    {
        for ( int idx=0; idx<nrknots; idx++ )
            ci += idx;
    }
    else
    {
        for( int idx=0; idx<nrknots; idx++ )
            ci += (nrknots-1-idx);
    }

    /*Triangulate: three consecutive vertices (idx0,idx,idx1) in polygon,
      remove cursize-2 Vertices, creating 1 triangle every time */
    int cursize = nrknots;
    int errcheck = 2*cursize;

    for( int idx=cursize-1; cursize>2; )
    {
        if ( 0 >= (errcheck--) )
            return false;

        const int idx0 = cursize<=idx ? 0 : idx;
        idx = cursize<=idx0+1 ? 0 : idx0+1;
        const int idx1 = cursize<=idx+1 ? 0 : idx+1;

        const Coord& pos0 = knots[ci[idx0]];
        const Coord& pos = knots[ci[idx]];
        const Coord& pos1 = knots[ci[idx1]];
        if ( (((pos.x-pos0.x)*(pos1.y-pos0.y)) -
              ((pos.y-pos0.y)*(pos1.x-pos0.x)))<0 )
            continue;

        bool isvalid = true;
        for ( int idy=0; idy<cursize; idy++ )
        {
            if( (idy==idx0) || (idy==idx) || (idy==idx1) )
                continue;

            if ( pointInTriangle2D(knots[ci[idy]],pos0,pos,pos1,0.0) )
            {
                isvalid = false;
                break;
            }
        }

        if ( isvalid )
        {
          res += ci[idx0];
          res += ci[idx];
          res += ci[idx1];

          /* remove idx from remaining polygon */
          for( int i=idx, j=idx+1; j<cursize; i++, j++ )
              ci[i] = ci[j];

          cursize--;
          errcheck = 2*cursize;
        }
    }

    return true;
}

#endif