GeoTessProfileNPoint.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 PROFILENPOINT_OBJECT_H
#define PROFILENPOINT_OBJECT_H
 
// **** _SYSTEM INCLUDES_ ******************************************************
 
#include <iostream>
#include <string>
#include <fstream>
#include <sstream>
 
// use standard library objects
using namespace std;
 
// **** _LOCAL INCLUDES_ *******************************************************
 
#include "GeoTessUtils.h"
#include "GeoTessException.h"
#include "GeoTessData.h"
#include "GeoTessProfile.h"
#include "GeoTessProfileType.h"
#include "GeoTessInterpolatorType.h"
 
// **** _BEGIN GEOTESS NAMESPACE_ **********************************************
 
namespace geotess {
 
// **** _FORWARD REFERENCES_ ***************************************************
 
class GeoTessMetaData;
class IFStreamAscii;
class IFStreamBinary;
 
// **** _CLASS DEFINITION_ *****************************************************
 
class GEOTESS_EXP_IMP GeoTessProfileNPoint : virtual public GeoTessProfile
{
private:
 
    int                                 nNodes;
 
    float*                          radii;
 
    GeoTessData**                           data;
 
    mutable double**        y2;
 
    int*                                pointIndices;
 
    GeoTessProfileNPoint() : GeoTessProfile(), nNodes(0), radii(NULL), data(NULL),
            y2(NULL), pointIndices(NULL) {};
 
    double*                         spline(float* x, GeoTessData** y, int attributeIndex,
            double yp1, double ypn) const;
 
    void                            check(int attributeIndex) const;
 
public:
 
    GeoTessProfileNPoint(float* r, GeoTessData** dat, int size) : GeoTessProfile(), nNodes(size), radii(NULL),
    data(NULL), y2(NULL), pointIndices(NULL)
    {
        radii = new float[size];
        data = new GeoTessData*[size];
        for (int i=0; i<size; ++i)
        {
            radii[i] = r[i];
            data[i] = dat[i];
        }
 
        if (radii[0] > radii[size - 1])
        {
            ostringstream os;
            os << endl << "ERROR in ProfileNPoint::ProfileNPoint" << endl
                    << "Profile has negative thickness" << endl;
            os << "radii = ";
            for (int i=0; i<size; ++i)
                os << radii[i] << ", ";
            os << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 4301);
        }
    }
 
    GeoTessProfileNPoint(const vector<float>& r, vector<GeoTessData*>& d)
    : GeoTessProfile(), nNodes(r.size()), radii(NULL),
    data(NULL), y2(NULL), pointIndices(NULL)
    {
        radii = new float[nNodes];
        data = new GeoTessData*[nNodes];
        for (int i=0; i<nNodes; ++i)
        {
            radii[i] = r[i];
            data[i] = d[i];
        }
 
        if (r.size() != d.size())
        {
            ostringstream os;
            os << endl << "ERROR in ProfileNPoint::ProfileNPoint" << endl
                    << "radii.size() != data.size()" << endl;
            os << "radii.size = " << r.size() << endl;
            os << "data.size   = " << d.size() << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 4302);
        }
 
        if (radii[0] > radii[nNodes - 1])
        {
            ostringstream os;
            os << endl << "ERROR in ProfileNPoint::ProfileNPoint" << endl
                    << "Profile has negative thickness" << endl;
            os << "radii = ";
            for (int i=0; i<nNodes; ++i)
                os << radii[i] << ", ";
            os << endl;
            throw GeoTessException(os, __FILE__, __LINE__, 4303);
        }
    }
 
    GeoTessProfileNPoint(float* rad, const vector<GeoTessData*>& dat);
 
    static  string          class_name() { return "ProfileNPoint"; };
 
    virtual int                 class_size() const
    { return (int) sizeof(GeoTessProfileNPoint); };
 
    virtual LONG_INT getMemory()
    {
        LONG_INT sz = (LONG_INT)sizeof(GeoTessProfileNPoint);
        if (nNodes > 0)
        {
            // add memory for radii
            sz += nNodes*(LONG_INT)sizeof(float);
 
            // add memory for data.  assumption is made that
            // GeoTessData objects at all nodes have the same type and
            // number of attributes.
            sz += nNodes*(LONG_INT)sizeof(GeoTessData*);
            sz += nNodes*data[0]->getMemory();
 
            if (pointIndices)
                sz += nNodes * (LONG_INT)(sizeof(int));
            if (y2)
                sz += data[0]->size() * nNodes * (LONG_INT)sizeof(double);
        }
        return sz;
    }
 
    virtual const GeoTessProfileType&   getType() const
    { return GeoTessProfileType::NPOINT; };
 
    virtual bool                operator == (const GeoTessProfile& p) const
                                                                                {
        if (!GeoTessProfile::operator==(p)) return false;
 
        if (nNodes != p.getNRadii()) return false;
 
        for (int i = 0; i < nNodes; ++i)
            if ((radii[i] != p.getRadius(i)) ||
                    (!(*(data[i]) == p.getData(i))))
                return false;
 
        return true;
                                                                                }
 
    bool                                isNaN(int nodeIndex, int attributeIndex)
    { return data[nodeIndex]->isNaN(attributeIndex); };
 
    virtual float               getRadius(int i) const { return radii[i]; };
 
    virtual int                 getNRadii() const { return nNodes; };
 
    virtual int                 getNData() const { return nNodes; };
 
    virtual float*          getRadii()
    {
        float* fa = new float [nNodes];
        for (int i=0; i<nNodes; ++i) fa[i] = radii[i];
        return fa;
    }
 
    virtual GeoTessData**           getData()
    {
        GeoTessData** da = new GeoTessData*[nNodes];
        for (int i=0; i<nNodes; ++i) da[i] = data[i];
        return da;
    }
 
 
    virtual GeoTessData*                getData(int i) { return data[i]; };
 
    virtual const GeoTessData&  getData(int i) const { return *(data[i]); };
 
    virtual void                setData(const vector<GeoTessData*>& inData);
 
    virtual void                    setRadii(const vector<float>& newRadii)
    { for (int i=0; i<nNodes; ++i) radii[i] = newRadii[i]; }
 
    virtual void setRadius(int index, float radius)
    { if (index < nNodes) radii[index] = radius; }
 
    virtual void                setData(int index, GeoTessData* inData)
    { delete data[index]; data[index] = inData; }
 
    virtual float               getRadiusTop() const { return radii[nNodes - 1]; };
 
    virtual const GeoTessData&  getDataTop() const { return *data[nNodes - 1]; };
 
    virtual GeoTessData*                getDataTop() { return data[nNodes - 1]; };
 
    virtual float               getRadiusBottom() const { return radii[0]; };
 
    virtual const GeoTessData&  getDataBottom() const { return *data[0]; };
 
    virtual GeoTessData*                getDataBottom() { return data[0]; };
 
 
    virtual double getValue(int attributeIndex, int radiusIndex) const
    {   return data[radiusIndex]->getDouble(attributeIndex); }
 
    virtual double getValueTop(int attributeIndex) const
    {   return data[nNodes-1]->getDouble(attributeIndex); }
 
    virtual double          getInterpolationCoefficient(int index, double radius) const;
 
    virtual double          getValue(const GeoTessInterpolatorType& radialType,
            int attributeIndex, double radius,
            bool allowRadiusOutOfRange) const;
 
    virtual int                 getRadiusIndex(double radius, int jlo) const;
 
    virtual double          getInterpolationCoefficient(int index, double radius,
            bool allowOutOfRange) const;
 
 
    GeoTessProfileNPoint(IFStreamBinary& ifs, GeoTessMetaData& gtmd);
 
    GeoTessProfileNPoint(IFStreamAscii& ifs, GeoTessMetaData& gtmd);
 
    virtual                         ~GeoTessProfileNPoint();
 
    virtual void                write(IFStreamBinary& ofs);
 
    virtual void                write(IFStreamAscii& ofs);
 
    virtual int                 findClosestRadiusIndex(double radius) const
    {
        int i = GeoTessProfile::getRadiusIndex(radius);
        return abs(radii[i+1] - radius) < abs(radii[i] - radius) ? i+1 : i;
    }
 
    virtual void                setPointIndex(int nodeIndex, int pointIndex)
    {
        if (pointIndices == NULL)
        {
            if (pointIndex < 0) return;
 
            pointIndices = new int [nNodes];
            for (int i = 0; i < nNodes; ++i) pointIndices[i] = -1;
        }
        pointIndices[nodeIndex] = pointIndex;
    }
 
    virtual void resetPointIndices()
    {
        if (pointIndices != NULL)
            delete[] pointIndices;
        pointIndices = NULL;
    }
 
    virtual int                 getPointIndex(int nodeIndex) const
    {   return pointIndices == NULL ? -1 : pointIndices[nodeIndex]; }
 
    virtual void                getWeights(map<int, double>& weights,
            double dkm, double radius, double hcoefficient) const
    {
        int node = GeoTessProfile::getRadiusIndex(radius);
 
        //TODO:  need getInterpolationCoefficient to work for cubic spline interpolator.  It currently does not.
 
        // get coefficient at node influencing position radius
 
        int pointIndex = pointIndices == NULL ? -1 : pointIndices[node];
 
        double c = getInterpolationCoefficient(node, radius, true);
        if (c > 0.0)
        {
            map<int, double>::iterator it = weights.find(pointIndex);
            if (it == weights.end())
                weights[pointIndex] = dkm * hcoefficient * c;
            else
                it->second += dkm * hcoefficient * c;
        }
 
        // now get coefficient at node+1  influencing position radius
 
        c = 1.0 - c;
        pointIndex = pointIndices == NULL ? -1 : pointIndices[node+1];
        if (c > 0.0)
        {
            map<int, double>::iterator it = weights.find(pointIndex);
            if (it == weights.end())
                weights[pointIndex] = dkm * hcoefficient * c;
            else
                it->second += dkm * hcoefficient * c;
        }
    }
 
    virtual void                    getCoefficients(map<int, double>& coefficients, double radius,
            double horizontalCoefficient) const
    {
        int node = GeoTessProfile::getRadiusIndex(radius);
 
        //TODO:  need getInterpolationCoefficient to work for cubic spline interpolator.  It currently does not.
        double c = getInterpolationCoefficient(node, radius, true);
 
        int pointIndex = pointIndices == NULL ? -1 : pointIndices[node];
 
        if (c > 0.0)
            coefficients[pointIndex] = c * horizontalCoefficient;
 
        pointIndex = pointIndices == NULL ? -1 : pointIndices[node+1];
 
        if (c < 1.0)
            coefficients[pointIndex] = (1.0 - c) * horizontalCoefficient;
    }
 
    virtual void setInterpolationCoefficients(const GeoTessInterpolatorType& interpType,
            vector<int>& nodeIndexes, vector<double>& coefficients,
            double& radius, bool& allowOutOfRange)
    {
        //TODO:  need this method to work for cubic spline interpolator.  It currently does not.
 
        if (radius < radii[0])
        {
            nodeIndexes.push_back(0);
            coefficients.push_back(allowOutOfRange ? 1 : NaN_DOUBLE);
        }
        else if (radius > radii[nNodes-1])
        {
            nodeIndexes.push_back(nNodes-1);
            coefficients.push_back(allowOutOfRange ? 1 : NaN_DOUBLE);
        }
        else
        {
            int index = getRadiusIndex(radius, -1);
            double c = ((double)radii[index + 1] - radius) /
                    ((double)radii[index + 1] - (double)radii[index]);
            nodeIndexes.push_back(index);
            coefficients.push_back(c);
            if (c < 1.)
            {
                nodeIndexes.push_back(index+1);
                coefficients.push_back(1.-c);
            }
        }
    }
 
    virtual GeoTessProfile*      copy()
    {
        GeoTessData** d = new GeoTessData* [nNodes];
        float* r = new float [nNodes];
        for (int i = 0; i < nNodes; ++i)
        {
            d[i] = data[i]->copy();
            r[i] = radii[i];
        }
        return new GeoTessProfileNPoint(r, d, nNodes);
    }
 
 
}; // end class ProfileNPoint
 
inline double   GeoTessProfileNPoint::getInterpolationCoefficient(int index,  double radius) const
{
    if (radius <= radii[index]) return 1.0;
    if (radius >= radii[index + 1]) return 0.0;
    return (radii[index + 1] - radius) / (radii[index + 1] - radii[index]);
}
 
inline double GeoTessProfileNPoint::getValue(const GeoTessInterpolatorType& radialType,
        int attributeIndex, double radius, bool allowRadiusOutOfRange) const
{
    if (!allowRadiusOutOfRange &&
            ((radius < (double)radii[0]) || (radius > (double)radii[nNodes-1])))
        return NaN_DOUBLE;
 
    int index = getRadiusIndex(radius, -1);
 
    double r0 = radii[index];
    double v0 = data[index]->getDouble(attributeIndex);
 
    if (radius <= r0) return v0;
 
    double r1 = radii[index + 1];
    double v1 = data[index + 1]->getDouble(attributeIndex);
 
    if (radius >= r1) return v1;
 
    double a = (r1 - radius) / (r1 - r0);
 
    double b = 1. - a;
    double v = a * v0 + b * v1;
 
    switch (radialType.ordinal())
    {
    case 0: // LINEAR
        return  v;
 
    case 2: // CUBIC_SPLINE
        check(attributeIndex);
        return v + ((a * a * a - a) * y2[attributeIndex][index]
                                                         + (b * b * b - b) * y2[attributeIndex][index + 1])
                                                         * (r1 - r0)    * (r1 - r0) / 6.0;
 
    default:
        ostringstream os;
        os << endl << "ERROR in ProfileNPoint::getValue" << endl
                << "InterpolatorType: " << radialType.name()
                << " cannot be applied to a Profile." << endl
                << "Must specify LINEAR or SPLINE" << endl;
        throw GeoTessException(os, __FILE__, __LINE__, 4304);
    }
}
 
 
inline void GeoTessProfileNPoint::check(int attributeIndex) const
{
    if (y2 == NULL)
    {
        y2 = new double* [data[0]->size()];
        for (int i=0; i<data[0]->size(); ++i)
            y2[i] = NULL;
    }
    if (y2[attributeIndex] == NULL)
        y2[attributeIndex] = spline(radii, data, attributeIndex, 1.0e30, 1.0e30);
}
 
 
} // end namespace geotess
 
#endif  // PROFILENPOINT_OBJECT_H