GeoTessPolygon.h

Go to the documentation of this file.

//- ****************************************************************************
//- 
//- Copyright 2009 Sandia Corporation. Under the terms of Contract
//- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
//- retains certain rights in this software.
//- 
//- BSD Open Source License.
//- All rights reserved.
//- 
//- Redistribution and use in source and binary forms, with or without
//- modification, are permitted provided that the following conditions are met:
//- 
//-    * Redistributions of source code must retain the above copyright notice,
//-      this list of conditions and the following disclaimer.
//-    * Redistributions in binary form must reproduce the above copyright
//-      notice, this list of conditions and the following disclaimer in the
//-      documentation and/or other materials provided with the distribution.
//-    * Neither the name of Sandia National Laboratories nor the names of its
//-      contributors may be used to endorse or promote products derived from
//-      this software without specific prior written permission.
//- 
//- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
//- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
//- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
//- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
//- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
//- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
//- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
//- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
//- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
//- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
//- POSSIBILITY OF SUCH DAMAGE.
//-
//- ****************************************************************************
 
#ifndef POLYGON_H_
#define POLYGON_H_
 
// **** _SYSTEM INCLUDES_ ******************************************************
 
#include <cstdio>
#include <iostream>
#include <iomanip>
#include <vector>
 
// use standard library objects
using namespace std;
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "CPPUtils.h"
#include "GeoTessUtils.h"
#include "GeoTessGreatCircle.h"
#include "IFStreamAscii.h"
 
// **** _BEGIN GEOTESS NAMESPACE_ **********************************************
 
namespace geotess {
 
// **** _FORWARD REFERENCES_ ***************************************************
 
// **** _CLASS DEFINITION_ *****************************************************
 
class GEOTESS_EXP_IMP GeoTessPolygon
{
protected:
 
    vector<GeoTessGreatCircle*> edges;
 
    double* referencePoint;
 
    bool referenceIn;
 
     static double TOLERANCE;
 
     bool global;
 
     bool lonFirst;
 
    void setup(vector<double*>& points);
 
    int edgeCrossings(GeoTessGreatCircle& gcRef);
 
    int refCount;
 
  public:
 
    virtual ~GeoTessPolygon();
 
    void* attachment;
 
    GeoTessPolygon();
 
    GeoTessPolygon(vector<double*>& points);
 
    GeoTessPolygon(const double* center, double radius, int nEdges);
 
    GeoTessPolygon(string filename);
 
    virtual string class_name() { return "Polygon"; };
 
    void addReference() { ++refCount; }
 
    void removeReference()
    {
        if (isNotReferenced())
        {
            ostringstream os;
            os << endl << "ERROR in Polygon::removeReference" << endl
                    << "Reference count (" << refCount << ") is already zero." << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 10001);
        }
 
        --refCount;
    }
 
    bool isNotReferenced() { return refCount == 0; }
 
    int getRefCount() { return refCount; }
 
    int size()
    {
        return edges.size();
    }
 
    static double getTolerance()
    {
        return TOLERANCE;
    }
 
    const double* const getReferencePoint()
    {
        return referencePoint;
    }
 
    const void getReferencePoint(double *u)
    {
        u[0] = referencePoint[0];
        u[1] = referencePoint[1];
        u[2] = referencePoint[2];
    }
 
    bool getReferencePointIn()
    {
        return referenceIn;
    }
 
    void invert()
    {
        referenceIn = !referenceIn;
    }
 
    void setReferencePoint(const double* refPoint, const bool &inside)
    {
        if (referencePoint == NULL)
            referencePoint = new double[3];
 
        referencePoint[0] = refPoint[0];
        referencePoint[1] = refPoint[1];
        referencePoint[2] = refPoint[2];
        referenceIn = inside;
    }
 
    void setReferencePoint(double lat, double lon, bool inside)
    {
        double r[3];
        GeoTessUtils::getVectorDegrees(lat, lon, r);
        setReferencePoint(r, inside);
    }
 
    bool containsAny(const vector<double*>& points)
    {
        for (int i=0; i<(int)points.size(); ++i)
            if (contains(points[i]))
                return true;
        return false;
    }
 
    bool containsAll(const vector<double*>& positions)
    {
        for (int i=0; i<(int)positions.size(); ++i)
            if (!contains(positions[i]))
                return false;
        return true;
    }
 
    bool contains(const double* x)
    {
        if (global || GeoTessUtils::dot(referencePoint, x) > cos(TOLERANCE))
            return referenceIn;
 
        GeoTessGreatCircle gcRef(referencePoint, x);
 
        return onBoundary(gcRef) || ((edgeCrossings(gcRef) % 2 == 0) == referenceIn);
    }
 
    bool onBoundary(GeoTessGreatCircle& gcRef)
    {
        // point being evaluated is gcRef.getLast()
        for (int i=0; i<(int)edges.size(); ++i)
        {
            // if point is very close to a polygon point, return true
            if (GeoTessUtils::dot(gcRef.getLast(), edges[i]->getFirst()) >= cos(TOLERANCE))
                return true;
 
            // if gcRef and edge are coincident, i.e., their normals are parallel or
            // anti-parallel, and the point being evaluated is between first and last
            // point of the edge, then return true;
            if (abs(GeoTessUtils::dot(gcRef.getNormal(), edges[i]->getNormal())) >= cos(TOLERANCE)
                    && edges[i]->getDistance(gcRef.getLast()) <= edges[i]->getDistance())
                return true;
        }
        return false;
    }
 
    bool onBoundary(const double* x)
    {
        GeoTessGreatCircle gc(referencePoint, x);
        return onBoundary(gc);
    }
 
    void getPoints(vector<double*> &points, const bool &repeatFirstPoint)
    {
        points.reserve(points.size()+edges.size()+1);
 
        for (int i = 0; i < (int)edges.size(); ++i)
        {
            double* point = new double[3];
            point[0] = edges[i]->getFirst()[0];
            point[1] = edges[i]->getFirst()[1];
            point[2] = edges[i]->getFirst()[2];
            points.push_back(point);
        }
 
        if (repeatFirstPoint)
        {
            double* point = new double[3];
            point[0] = edges[0]->getFirst()[0];
            point[1] = edges[0]->getFirst()[1];
            point[2] = edges[0]->getFirst()[2];
            points.push_back(point);
        }
    }
 
    void getPoints(vector<double*> &points, const bool &repeatFirstPoint, const double &maxSpacing)
    {
        points.reserve(points.size()+edges.size()+1);
 
        int n;
        double dx;
        for (int i = 0; i < (int)edges.size(); ++i)
        {
            n = (int)ceil(edges[i]->getDistance()/maxSpacing);
            dx = edges[i]->getDistance()/n;
            for (int j=0; j<n; ++j)
                points.push_back(edges[i]->getPoint(j*dx));
        }
 
        if (repeatFirstPoint)
        {
            double* point = new double[3];
            point[0] = edges[0]->getFirst()[0];
            point[1] = edges[0]->getFirst()[1];
            point[2] = edges[0]->getFirst()[2];
            points.push_back(point);
        }
    }
 
    const double* getPoint(int index)
    {
        return edges[index]->getFirst();
    }
 
    double getArea()
    {
        double a, area = 0;
        GeoTessGreatCircle* edge;
        GeoTessGreatCircle* previous = edges[edges.size()-1];
 
        for (int i=0; i<(int)edges.size(); ++i)
        {
            edge = edges[i];
            a = PI - GeoTessUtils::angle(previous->getNormal(), edge->getNormal());
 
            if (GeoTessUtils::scalarTripleProduct(previous->getNormal(), edge->getNormal(),
                    edge->getFirst()) < 0)
                area += a;
            else
                area += 2*PI - a;
 
            previous = edge;
        }
        area -= (edges.size()-2)*PI;
        return area;
    }
 
    double getAreaSmall()
    {
        double area = getArea();
        return area <= 2*PI ? area : 4*PI-area;
    }
 
    double getAreaLarge()
    {
        double area = getArea();
        return area >= 2*PI ? area : 4*PI-area;
    }
 
    string str(const bool& repeatFirstPoint, const bool& latFirst,
            const double& minLongitude = -180);
 
    virtual void write(const string& outputFileName);
 
 
    virtual void loadAscii(vector<string>& records);
 
 
}; // end class Polygon
 
} // end namespace geotess
 
#endif /* POLYGON_H_ */