GeoTessPosition.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 GEOTESSPOSITION_OBJECT_H
#define GEOTESSPOSITION_OBJECT_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
 
#include <iostream>
#include <string>
#include <fstream>
#include <vector>
#include <sstream>
 
// use standard library objects
using namespace std;
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "CPPUtils.h"
#include "GeoTessUtils.h"
#include "GeoTessModel.h"
#include "GeoTessGrid.h"
#include "GeoTessMetaData.h"
#include "CpuTimer.h"
#include "GeoTessProfile.h"
#include "GeoTessData.h"
 
// **** _BEGIN GEOTESS NAMESPACE_ **********************************************
 
namespace geotess {
 
// **** _FORWARD REFERENCES_ ***************************************************
 
class GeoTessInterpolatorType;
class GeoTessProfile;
 
// **** _CLASS DEFINITION_ *****************************************************
 
class GEOTESS_EXP_IMP GeoTessPosition
{
private:
 
    int refCount;
 
    int*                                    maxTessLevel;
 
    int*                                    tessLevels;
 
    int*                                    triangle;
 
    double                              errorValue;
 
    void                                    updatePosition2D(int layerid, const double* const uVector);
 
    void                                    checkTessellation(const int& tid)
    {
        if (triangle[tid] < 0)
        {
            tessLevels[tid] = 0;
            triangle[tid] = grid.getTriangle(tid, 0, 0);
            getContainingTriangle(tid);
            update2D(tid);
        }
    }
 
    void    updateRadius(int lid, double rad)
    {
        if ((radius < 0.) || (rad != radius) || (lid != layerId))
        {
            radius = rad;
            layerId = lid;
            tessid = layerTessIds[lid];
            checkTessellation(tessid);
            clearRadialCoefficients();
        }
    }
 
    void clearRadialCoefficients()
    {
        for (int i=0; i<(int)radialIndexes.size(); ++i)
        {
            radialIndexes[i].clear();
            radialCoefficients[i].clear();
        }
    }
 
    void updateRadialCoefficients()
    {
        // make sure dimensions of radialIndexes and radialCoefficients
        // are at least as big as the number of vertices involved in interpolation.
        if (radialIndexes.size() < vertices[tessid].size())
        {
            radialIndexes.resize(vertices[tessid].size());
            radialCoefficients.resize(vertices[tessid].size());
        }
 
        // if radial interpolation coefficients have already been computed
        // do nothing.
        if (radialIndexes[0].size() == 0)
        {
            vector<int>& v = vertices[tessid];
            for (int i = 0; i < (int)v.size(); ++i)
                modlProfiles[v[i]][layerId]->setInterpolationCoefficients(radialInterpolatorType,
                        radialIndexes[i], radialCoefficients[i], radius, radiusOutOfRangeAllowed);
        }
    }
 
    static const double TWALK_TOLERANCE;
 
    void                                    getContainingTriangle(int tid);
 
    int biggerTriangle(int t1, int t2)
    {
        const int* tv = grid.getTriangleVertexIndexes(t1);
        double dot1 = GeoTessUtils::dot(grid.getVertex(tv[0]), grid.getVertex(tv[1]))
                + GeoTessUtils::dot(grid.getVertex(tv[1]), grid.getVertex(tv[2]))
                + GeoTessUtils::dot(grid.getVertex(tv[2]), grid.getVertex(tv[0]));
 
        tv = grid.getTriangleVertexIndexes(t2);
        double dot2 = GeoTessUtils::dot(grid.getVertex(tv[0]), grid.getVertex(tv[1]))
                + GeoTessUtils::dot(grid.getVertex(tv[1]), grid.getVertex(tv[2]))
                + GeoTessUtils::dot(grid.getVertex(tv[2]), grid.getVertex(tv[0]));
 
        return dot2 > dot1 ? t1 : t2;
 
    }
 
protected:
 
 
    double                              radius;
 
    double                              earthRadius;
 
    vector<double>              layerRadii;
 
    int                                     layerId;
 
    int                                     tessid;
 
    vector< vector< int> >  vertices;
 
    vector< vector< double> >   linearCoefficients;
 
    //vector< vector< double> >     nnCoefficients;
 
    vector< vector< double> >   hCoefficients;
 
    vector<vector<int> > radialIndexes;
 
    vector<vector<double> > radialCoefficients;
 
 
    const GeoTessInterpolatorType& radialInterpolatorType;
 
    GeoTessModel*       model;
 
    GeoTessGrid&        grid;
 
    GeoTessProfile ***  modlProfiles;
    double const* const*          gridVertices;
    int const* const*             gridTriangles;
    int const*                    gridDescendants;
    const vector<vector<Edge*> >& gridEdges;
    const int*                    layerTessIds;
    int                           nLayers;
 
    bool                                                radiusOutOfRangeAllowed;
 
    double                              unitVector[3];
 
    int                                     getTriangle(int tid)
    { checkTessellation(tid); return triangle[tid]; };
 
    virtual void                    update2D(int tid) = ABSTRACT;
 
    GeoTessPosition(GeoTessModel* model, const GeoTessInterpolatorType& radialType);
 
    long memory()
    {
        long m = 0;
 
        // memory for int* triangle, tessLevels and maxTessLevel.
        m += (long)(grid.getNTessellations() * 3 * sizeof(int));
 
        // vector<double> layerRadii;
        m += (long)(layerRadii.capacity() * sizeof(double));
 
        // vector< vector< int> > vertices;
        m += (long)(vertices.capacity()*sizeof(vector<int>));
        for (int i=0; i<(int)vertices.size(); ++i)
            m += (long)(vertices[i].capacity()*sizeof(int));
 
        // vector< vector< double> > linearCoefficients;
        m += (long)(linearCoefficients.capacity()*sizeof(vector<double>));
        for (int i=0; i<(int)linearCoefficients.size(); ++i)
            m += (long)(linearCoefficients[i].capacity()*sizeof(double));
 
        // vector< vector< double> > hCoefficients;
        m += (long)(hCoefficients.capacity()*sizeof(vector<double>));
        for (int i=0; i<(int)hCoefficients.size(); ++i)
            m += (long)(hCoefficients[i].capacity()*sizeof(double));
 
        // vector<vector<int> > radialIndexes;
        m += (long)(radialIndexes.capacity()*sizeof(vector<int>));
        for (int i=0; i<(int)radialIndexes.size(); ++i)
            m += (long)(radialIndexes[i].capacity()*sizeof(int));
 
        // vector<vector<double> > radialCoefficients;
        m += (long)(radialCoefficients.capacity()*sizeof(vector<double>));
        for (int i=0; i<(int)radialCoefficients.size(); ++i)
            m += (long)(radialCoefficients[i].capacity()*sizeof(double));
 
        return m;
    }
 
 
public:
 
    static GeoTessPosition* getGeoTessPosition(GeoTessModel* model);
 
    static GeoTessPosition* getGeoTessPosition(GeoTessModel* model,
            const GeoTessInterpolatorType& horizontalType);
 
    static GeoTessPosition* getGeoTessPosition(GeoTessModel* model,
            const GeoTessInterpolatorType& horizontalType, const GeoTessInterpolatorType& radialType);
 
    virtual                             ~GeoTessPosition();
 
    virtual const GeoTessInterpolatorType&  getInterpolatorType() const = ABSTRACT;
 
    virtual long getMemory() = ABSTRACT;
 
    virtual double              getValue(int attribute);
 
    void                                    setModel(GeoTessModel* newModel);
 
    bool                                    noEmptyProfiles();
 
    void                                    set(double lat, double lon, double depth)
    {
        double uVector[3] = {0.0, 0.0, 0.0};
        GeoTessUtils::getVectorDegrees(lat, lon, uVector);
        double newRadius = GeoTessUtils::getEarthRadius(uVector) - depth;
        set(uVector, newRadius);
    }
 
    void set(const double* const uVector, const double& newRadius)
    {
        updatePosition2D(nLayers - 1, uVector);
 
        int lid = getLayerId(newRadius);
 
        updatePosition2D(lid, uVector);
        updateRadius(lid, newRadius);
    }
 
    void set(int layid, double lat, double lon, double depth)
    {
        if (layid < 0)
            set(lat, lon, depth);
        else
        {
            double uVector[3];
            GeoTessUtils::getVectorDegrees(lat, lon, uVector);
            updatePosition2D(layid, uVector);
            updateRadius(layid, GeoTessUtils::getEarthRadius(uVector) - depth);
        }
    }
 
    void                                    set(int layid, const double* const uVector, double rad)
    {
        if (layid < 0)
            set(uVector, rad);
        else
        {
            updatePosition2D(layid, uVector);
            if (rad >= 0.0) updateRadius(layid, rad);
        }
    }
 
    void                                    setTop(int layid, const double* const uVector)
    {
        updatePosition2D(layid, uVector);
        updateRadius(layid, getRadiusTop(layid));
    }
 
    void                                    setBottom(int layid, const double* const uVector)
    {
        updatePosition2D(layid, uVector);
        updateRadius(layid, getRadiusBottom(layid));
    }
 
    void                                    setRadius(int layid, double rad)
    {
        if (tessid < 0)
        {
            ostringstream os;
            os << endl << "ERROR in GeoTessPosition::setRadius" << endl
                    << "Geographic position has not been specified." << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 3001);
        }
        updateRadius(layid, rad);
    }
 
    void setRadius(double rad)
    {
        if (tessid < 0)
            if (tessid < 0)
            {
                ostringstream os;
                os << endl << "ERROR in GeoTessPosition::setRadius" << endl
                        << "Geographic position has not been specified." << endl;
                throw GeoTessException(os, __FILE__, __LINE__, 3002);
            }
        updateRadius(getLayerId(rad), rad);
    }
 
    void setDepth(int layer, double depth)
    {
        setRadius(layer, getEarthRadius()-depth);
    }
 
    void setDepth(double depth)
    {
        setRadius(getEarthRadius()-depth);
    }
 
    void                                    setTop(int layid)
    {
        if (tessid < 0)
        {
            ostringstream os;
            os << endl << "ERROR in GeoTessPosition::setRadius" << endl
                    << "Geographic position has not been specified." << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 3003);
        }
        tessid = layerTessIds[layid];
        checkTessellation(tessid);
        updateRadius(layid, getRadiusTop(layid));
    }
 
    void                                    setBottom(int layid)
    {
        if (tessid < 0)
        {
            ostringstream os;
            os << endl << "ERROR in GeoTessPosition::setRadius" << endl
                    << "Geographic position has not been specified." << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 3004);
        }
        tessid = layerTessIds[layid];
        checkTessellation(tessid);
        updateRadius(layid, getRadiusBottom(layid));
    }
 
    double                              getRadiusTop(int layid);
 
    double                              getRadiusBottom(int layid);
 
    double                              getEarthRadius()
    {
        if (earthRadius < 0.)
            earthRadius = GeoTessUtils::getEarthRadius(unitVector);
        return earthRadius;
    }
 
    double*                             getVector() { return unitVector; };
 
    void                                    copyVector(double* u)
    { u[0] = unitVector[0]; u[1] = unitVector[1]; u[2] = unitVector[2]; };
 
    int                                     getTriangle() { return triangle[tessid]; };
 
    int                                     getNVertices() { return vertices[tessid].size(); };
 
    const vector<int>& getVertices() const { return vertices[tessid]; };
 
    int getIndexOfClosestVertex() const;
 
    const double*         getClosestVertex() const
    {
        return grid.getVertex(getIndexOfClosestVertex());
    }
 
    int                                     getVertex(int index)
    {   return vertices[tessid][index]; }
 
    void    getCoefficients(map<int, double>& coefficients)
    {
        coefficients.clear();
        vector<int>& vtx = vertices[tessid];
        vector<double>& hc = hCoefficients[tessid];
        for (int i = 0; i < (int) vtx.size(); ++i)
            modlProfiles[vtx[i]][layerId]->getCoefficients(coefficients, radius, hc[i]);
    }
 
    void                                    setMaxTessLevel(int layid, int maxTess)
    {
        maxTessLevel[layerTessIds[layid]] = maxTess;
        triangle[layerTessIds[layid]] = -1;
        checkTessellation(tessid);
    }
 
    int                                     getMaxTessLevel(int layid)
    {   return maxTessLevel[layerTessIds[layid]]; }
 
    int                                     getTessLevel() const { return tessLevels[tessid]; };
 
    int getTessLevel(const int& tId) { checkTessellation(tId); return tessLevels[tId]; }
 
    double                              getRadiusTop()
    {   return getRadiusTop(layerId); };
 
    double                              getRadiusBottom()
    {   return getRadiusBottom(layerId); };
 
    double                              getDepthTop()
    {
        return  getEarthRadius() - getRadiusTop(layerId);
    }
 
    double                              getDepthBottom()
    {
        return  getEarthRadius() - getRadiusBottom(layerId);
    }
 
    double                              getDepthTop(int layid)
    {
        return  getEarthRadius() - getRadiusTop(layid);
    }
 
    double                              getDepthBottom(int layid)
    {
        return  getEarthRadius() - getRadiusBottom(layid);
    }
 
    double                              getLayerThickness(int layid)
    {
        return  getRadiusTop(layid) - getRadiusBottom(layid);
    }
 
    double                              getLayerThickness()
    {
        return  getRadiusTop() - getRadiusBottom();
    }
 
    double                              getRadius() { return radius; };
 
    double                              getDepth()
    {   return getEarthRadius() - radius; };
 
    GeoTessModel* getModel() { return model; }
 
    int getTessID() { return tessid; };
 
    int                                     getLayerId(double rad)
    {
        for (int i=0; i<nLayers; ++i)
            if (rad <= getRadiusTop(i))
                return i;
 
        for (int i=nLayers-1; i>=0; --i)
            if (getLayerThickness(i) > 0)
                return i;
 
        return nLayers-1;
    };
 
    int                                     getLayerId() { return layerId; };
 
    string                              toString();
 
    double                              getErrorValue()
    {
        return errorValue;
    }
 
    void                                    setErrorValue(double errVal)
    {
        errorValue = errVal;
    }
 
    int                                   getVertexIndex()
    {
        vector<int>& vtid = vertices[tessid];
        vector<double>& htid = hCoefficients[tessid];
        for (int v = 0; v < (int) vtid.size(); ++v)
            if (htid[v] > 0.999999999)
                return vtid[v];
        return -1;
    }
 
    void getWeights(map<int, double>& weights, double dkm)
    {
        vector<int>& vtx = vertices[tessid];
        vector<double>& htid = hCoefficients[tessid];
        model->getPointMap();
        for (int i = 0; i < (int) vtx.size(); ++i)
            modlProfiles[vtx[i]][layerId]->getWeights(weights, dkm, radius, htid[i]);
    }
 
    const vector<double>& getHorizontalCoefficients() const
    {
        return hCoefficients[tessid];
    }
 
    double getHorizontalCoefficient(int index) const
    {
        return hCoefficients[tessid][index];
    }
 
    bool isRadiusOutOfRangeAllowed()
    {
        return radiusOutOfRangeAllowed;
    }
 
    void setRadiusOutOfRangeAllowed(bool allowed)
    {
        if (allowed != radiusOutOfRangeAllowed)
            clearRadialCoefficients();
 
        radiusOutOfRangeAllowed = allowed;
    }
 
    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__, 3005);
        }
 
        --refCount;
    }
 
    int getReferenceCount()
    {
        return refCount;
    }
 
    bool isNotReferenced() { return refCount == 0; }
 
}; // end class GeoTessModel
 
} // end namespace geotess
 
#endif  // GEOTESSPOSITION_OBJECT_H