GeoTessCPP 2.6.1 Software to facilitate storage and retrieval of 3D information about the Earth.

Loading...
Searching...
No Matches
GeoTessProfileNPoint.h
Go to the documentation of this file.
1//- ****************************************************************************
2//-
3//- Copyright 2009 Sandia Corporation. Under the terms of Contract
4//- DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
5//- retains certain rights in this software.
6//-
7//- BSD Open Source License.
8//- All rights reserved.
9//-
10//- Redistribution and use in source and binary forms, with or without
11//- modification, are permitted provided that the following conditions are met:
12//-
13//- * Redistributions of source code must retain the above copyright notice,
14//- this list of conditions and the following disclaimer.
15//- * Redistributions in binary form must reproduce the above copyright
16//- notice, this list of conditions and the following disclaimer in the
17//- documentation and/or other materials provided with the distribution.
18//- * Neither the name of Sandia National Laboratories nor the names of its
19//- contributors may be used to endorse or promote products derived from
20//- this software without specific prior written permission.
21//-
22//- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
23//- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
24//- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
25//- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
26//- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
27//- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
28//- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
29//- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
30//- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
31//- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
32//- POSSIBILITY OF SUCH DAMAGE.
33//-
34//- ****************************************************************************
35
36#ifndef PROFILENPOINT_OBJECT_H
37#define PROFILENPOINT_OBJECT_H
38
39// **** _SYSTEM INCLUDES_ ******************************************************
40
41#include <iostream>
42#include <string>
43#include <fstream>
44#include <sstream>
45
46// use standard library objects
47 using namespace std;
48
49// **** _LOCAL INCLUDES_ *******************************************************
50
51#include "GeoTessUtils.h"
52#include "GeoTessException.h"
53#include "GeoTessData.h"
54#include "GeoTessProfile.h"
55#include "GeoTessProfileType.h"
56#include "GeoTessInterpolatorType.h"
57
58// **** _BEGIN GEOTESS NAMESPACE_ **********************************************
59
60 namespace geotess {
61
62// **** _FORWARD REFERENCES_ ***************************************************
63
64 class GeoTessMetaData;
65 class IFStreamAscii;
66 class IFStreamBinary;
67
68// **** _CLASS DEFINITION_ *****************************************************
69
80 class GEOTESS_EXP_IMP GeoTessProfileNPoint : virtual public GeoTessProfile
81{
82 private:
83
87 int nNodes;
88
92 float* radii;
93
97 GeoTessData** data;
98
105 mutable double** y2;
106
114 int* pointIndices;
115
119 GeoTessProfileNPoint() : GeoTessProfile(), nNodes(0), radii(NULL), data(NULL),
120 y2(NULL), pointIndices(NULL) {};
121
126 double* spline(float* x, GeoTessData** y, int attributeIndex,
127 double yp1, double ypn) const;
128
134 void check(int attributeIndex) const;
135
136 public:
137
150 GeoTessProfileNPoint(float* r, GeoTessData** dat, int size) : GeoTessProfile(), nNodes(size), radii(NULL),
151 data(NULL), y2(NULL), pointIndices(NULL)
152 {
153 radii = new float[size];
154 data = new GeoTessData*[size];
155 for (int i=0; i<size; ++i)
156 {
157 radii[i] = r[i];
158 data[i] = dat[i];
159 }
160
161 if (radii[0] > radii[size - 1])
162 {
163 ostringstream os;
164 os << endl << "ERROR in ProfileNPoint::ProfileNPoint" << endl
165 << "Profile has negative thickness" << endl;
166 os << "radii = ";
167 for (int i=0; i<size; ++i)
168 os << radii[i] << ", ";
169 os << endl;
170 throw GeoTessException(os, __FILE__, __LINE__, 4301);
171 }
172 }
173
185 GeoTessProfileNPoint(const vector<float>& r, vector<GeoTessData*>& d)
186 : GeoTessProfile(), nNodes(r.size()), radii(NULL),
187 data(NULL), y2(NULL), pointIndices(NULL)
188 {
189 radii = new float[nNodes];
190 data = new GeoTessData*[nNodes];
191 for (int i=0; i<nNodes; ++i)
192 {
193 radii[i] = r[i];
194 data[i] = d[i];
195 }
196
197 if (r.size() != d.size())
198 {
199 ostringstream os;
200 os << endl << "ERROR in ProfileNPoint::ProfileNPoint" << endl
201 << "radii.size() != data.size()" << endl;
202 os << "radii.size = " << r.size() << endl;
203 os << "data.size = " << d.size() << endl;
204 throw GeoTessException(os, __FILE__, __LINE__, 4302);
205 }
206
207 if (radii[0] > radii[nNodes - 1])
208 {
209 ostringstream os;
210 os << endl << "ERROR in ProfileNPoint::ProfileNPoint" << endl
211 << "Profile has negative thickness" << endl;
212 os << "radii = ";
213 for (int i=0; i<nNodes; ++i)
214 os << radii[i] << ", ";
215 os << endl;
216 throw GeoTessException(os, __FILE__, __LINE__, 4303);
217 }
218 }
219
225 GeoTessProfileNPoint(float* rad, const vector<GeoTessData*>& dat);
226
230 static string class_name() { return "ProfileNPoint"; };
231
235 virtual int class_size() const
236 { return (int) sizeof(GeoTessProfileNPoint); };
237
238 virtual LONG_INT getMemory()
239 {
240 LONG_INT sz = (LONG_INT)sizeof(GeoTessProfileNPoint);
241 if (nNodes > 0)
242 {
243 // add memory for radii
244 sz += nNodes*(LONG_INT)sizeof(float);
245
246 // add memory for data. assumption is made that
247 // GeoTessData objects at all nodes have the same type and
248 // number of attributes.
249 sz += nNodes*(LONG_INT)sizeof(GeoTessData*);
250 sz += nNodes*data[0]->getMemory();
251
252 if (pointIndices)
253 sz += nNodes * (LONG_INT)(sizeof(int));
254 if (y2)
255 sz += data[0]->size() * nNodes * (LONG_INT)sizeof(double);
256 }
257 return sz;
258 }
259
265 virtual const GeoTessProfileType& getType() const
266 { return GeoTessProfileType::NPOINT; };
267
271 virtual bool operator == (const GeoTessProfile& p) const
272 {
273 if (!GeoTessProfile::operator==(p)) return false;
274
275 if (nNodes != p.getNRadii()) return false;
276
277 for (int i = 0; i < nNodes; ++i)
278 if ((radii[i] != p.getRadius(i)) ||
279 (!(*(data[i]) == p.getData(i))))
280 return false;
281
282 return true;
283 }
284
292 bool isNaN(int nodeIndex, int attributeIndex)
293 { return data[nodeIndex]->isNaN(attributeIndex); };
294
299 virtual float getRadius(int i) const { return radii[i]; };
300
304 virtual int getNRadii() const { return nNodes; };
305
309 virtual int getNData() const { return nNodes; };
310
315 virtual float* getRadii()
316 {
317 float* fa = new float [nNodes];
318 for (int i=0; i<nNodes; ++i) fa[i] = radii[i];
319 return fa;
320 }
321
327 virtual GeoTessData** getData()
328 {
329 GeoTessData** da = new GeoTessData*[nNodes];
330 for (int i=0; i<nNodes; ++i) da[i] = data[i];
331 return da;
332 }
333
334
338 virtual GeoTessData* getData(int i) { return data[i]; };
339
343 virtual const GeoTessData& getData(int i) const { return *(data[i]); };
344
348 virtual void setData(const vector<GeoTessData*>& inData);
349
353 virtual void setRadii(const vector<float>& newRadii)
354 { for (int i=0; i<nNodes; ++i) radii[i] = newRadii[i]; }
355
356 virtual void setRadius(int index, float radius)
357 { if (index < nNodes) radii[index] = radius; }
358
362 virtual void setData(int index, GeoTessData* inData)
363 { delete data[index]; data[index] = inData; }
364
368 virtual float getRadiusTop() const { return radii[nNodes - 1]; };
369
373 virtual const GeoTessData& getDataTop() const { return *data[nNodes - 1]; };
374
378 virtual GeoTessData* getDataTop() { return data[nNodes - 1]; };
379
383 virtual float getRadiusBottom() const { return radii[0]; };
384
388 virtual const GeoTessData& getDataBottom() const { return *data[0]; };
389
393 virtual GeoTessData* getDataBottom() { return data[0]; };
394
395
400 virtual double getValue(int attributeIndex, int radiusIndex) const
401 { return data[radiusIndex]->getDouble(attributeIndex); }
402
409 virtual double getValueTop(int attributeIndex) const
410 { return data[nNodes-1]->getDouble(attributeIndex); }
411
416 virtual double getInterpolationCoefficient(int index, double radius) const;
417
422 virtual double getValue(const GeoTessInterpolatorType& radialType,
423 int attributeIndex, double radius,
424 bool allowRadiusOutOfRange) const;
425
430 virtual int getRadiusIndex(double radius, int jlo) const;
431
435 virtual double getInterpolationCoefficient(int index, double radius,
436 bool allowOutOfRange) const;
437
439
443 GeoTessProfileNPoint(IFStreamBinary& ifs, GeoTessMetaData& gtmd);
444
448 GeoTessProfileNPoint(IFStreamAscii& ifs, GeoTessMetaData& gtmd);
449
453 virtual ~GeoTessProfileNPoint();
454
458 virtual void write(IFStreamBinary& ofs);
459
463 virtual void write(IFStreamAscii& ofs);
464
473 virtual int findClosestRadiusIndex(double radius) const
474 {
475 int i = GeoTessProfile::getRadiusIndex(radius);
476 return abs(radii[i+1] - radius) < abs(radii[i] - radius) ? i+1 : i;
477 }
478
486 virtual void setPointIndex(int nodeIndex, int pointIndex)
487 {
488 if (pointIndices == NULL)
489 {
490 if (pointIndex < 0) return;
491
492 pointIndices = new int [nNodes];
493 for (int i = 0; i < nNodes; ++i) pointIndices[i] = -1;
494 }
495 pointIndices[nodeIndex] = pointIndex;
496 }
497
505 virtual void resetPointIndices()
506 {
507 if (pointIndices != NULL)
508 delete[] pointIndices;
509 pointIndices = NULL;
510 }
511
519 virtual int getPointIndex(int nodeIndex) const
520 { return pointIndices == NULL ? -1 : pointIndices[nodeIndex]; }
521
528 virtual void getWeights(map<int, double>& weights,
529 double dkm, double radius, double hcoefficient) const
530 {
531 int node = GeoTessProfile::getRadiusIndex(radius);
532
533 //TODO: need getInterpolationCoefficient to work for cubic spline interpolator. It currently does not.
534
535 // get coefficient at node influencing position radius
536
537 int pointIndex = pointIndices == NULL ? -1 : pointIndices[node];
538
539 double c = getInterpolationCoefficient(node, radius, true);
540 if (c > 0.0)
541 {
542 map<int, double>::iterator it = weights.find(pointIndex);
543 if (it == weights.end())
544 weights[pointIndex] = dkm * hcoefficient * c;
545 else
546 it->second += dkm * hcoefficient * c;
547 }
548
549 // now get coefficient at node+1 influencing position radius
550
551 c = 1.0 - c;
552 pointIndex = pointIndices == NULL ? -1 : pointIndices[node+1];
553 if (c > 0.0)
554 {
555 map<int, double>::iterator it = weights.find(pointIndex);
556 if (it == weights.end())
557 weights[pointIndex] = dkm * hcoefficient * c;
558 else
559 it->second += dkm * hcoefficient * c;
560 }
561 }
562
569 virtual void getCoefficients(map<int, double>& coefficients, double radius,
570 double horizontalCoefficient) const
571 {
572 int node = GeoTessProfile::getRadiusIndex(radius);
573
574 //TODO: need getInterpolationCoefficient to work for cubic spline interpolator. It currently does not.
575 double c = getInterpolationCoefficient(node, radius, true);
576
577 int pointIndex = pointIndices == NULL ? -1 : pointIndices[node];
578
579 if (c > 0.0)
580 coefficients[pointIndex] = c * horizontalCoefficient;
581
582 pointIndex = pointIndices == NULL ? -1 : pointIndices[node+1];
583
584 if (c < 1.0)
585 coefficients[pointIndex] = (1.0 - c) * horizontalCoefficient;
586 }
587
588 virtual void setInterpolationCoefficients(const GeoTessInterpolatorType& interpType,
589 vector<int>& nodeIndexes, vector<double>& coefficients,
590 double& radius, bool& allowOutOfRange)
591 {
592 //TODO: need this method to work for cubic spline interpolator. It currently does not.
593
594 if (radius < radii[0])
595 {
596 nodeIndexes.push_back(0);
597 coefficients.push_back(allowOutOfRange ? 1 : NaN_DOUBLE);
598 }
599 else if (radius > radii[nNodes-1])
600 {
601 nodeIndexes.push_back(nNodes-1);
602 coefficients.push_back(allowOutOfRange ? 1 : NaN_DOUBLE);
603 }
604 else
605 {
606 int index = getRadiusIndex(radius, -1);
607 double c = ((double)radii[index + 1] - radius) /
608 ((double)radii[index + 1] - (double)radii[index]);
609 nodeIndexes.push_back(index);
610 coefficients.push_back(c);
611 if (c < 1.)
612 {
613 nodeIndexes.push_back(index+1);
614 coefficients.push_back(1.-c);
615 }
616 }
617 }
618
622 virtual GeoTessProfile* copy()
623 {
624 GeoTessData** d = new GeoTessData* [nNodes];
625 float* r = new float [nNodes];
626 for (int i = 0; i < nNodes; ++i)
627 {
628 d[i] = data[i]->copy();
629 r[i] = radii[i];
630 }
631 return new GeoTessProfileNPoint(r, d, nNodes);
632 }
633
635
636}; // end class ProfileNPoint
637
649 inline double GeoTessProfileNPoint::getInterpolationCoefficient(int index, double radius) const
650{
651 if (radius <= radii[index]) return 1.0;
652 if (radius >= radii[index + 1]) return 0.0;
653 return (radii[index + 1] - radius) / (radii[index + 1] - radii[index]);
654}
655
668 inline double GeoTessProfileNPoint::getValue(const GeoTessInterpolatorType& radialType,
669 int attributeIndex, double radius, bool allowRadiusOutOfRange) const
670{
671 if (!allowRadiusOutOfRange &&
672 ((radius < (double)radii[0]) || (radius > (double)radii[nNodes-1])))
673 return NaN_DOUBLE;
674
675 int index = getRadiusIndex(radius, -1);
676
677 double r0 = radii[index];
678 double v0 = data[index]->getDouble(attributeIndex);
679
680 if (radius <= r0) return v0;
681
682 double r1 = radii[index + 1];
683 double v1 = data[index + 1]->getDouble(attributeIndex);
684
685 if (radius >= r1) return v1;
686
687 double a = (r1 - radius) / (r1 - r0);
688
689 double b = 1. - a;
690 double v = a * v0 + b * v1;
691
692 switch (radialType.ordinal())
693 {
694 case 0: // LINEAR
695 return v;
696
697 case 2: // CUBIC_SPLINE
698 check(attributeIndex);
699 return v + ((a * a * a - a) * y2[attributeIndex][index]
700 + (b * b * b - b) * y2[attributeIndex][index + 1])
701 * (r1 - r0) * (r1 - r0) / 6.0;
702
703 default:
704 ostringstream os;
705 os << endl << "ERROR in ProfileNPoint::getValue" << endl
706 << "InterpolatorType: " << radialType.name()
707 << " cannot be applied to a Profile." << endl
708 << "Must specify LINEAR or SPLINE" << endl;
709 throw GeoTessException(os, __FILE__, __LINE__, 4304);
710 }
711}
712
714
720 inline void GeoTessProfileNPoint::check(int attributeIndex) const
721{
722 if (y2 == NULL)
723 {
724 y2 = new double* [data[0]->size()];
725 for (int i=0; i<data[0]->size(); ++i)
726 y2[i] = NULL;
727 }
728 if (y2[attributeIndex] == NULL)
729 y2[attributeIndex] = spline(radii, data, attributeIndex, 1.0e30, 1.0e30);
730}
731
733
734} // end namespace geotess
735
736#endif // PROFILENPOINT_OBJECT_H
GEOTESS_EXP_IMP
#define GEOTESS_EXP_IMP
Definition CPPGlobals.h:71
LONG_INT
#define LONG_INT
Definition CPPGlobals.h:111
GeoTessData.h
GeoTessException.h
GeoTessInterpolatorType.h
GeoTessProfile.h
GeoTessProfileType.h
GeoTessUtils.h
geotess::GeoTessData
Abstract base class that manages the data values attached to a single grid point.
Definition GeoTessData.h:76
geotess::GeoTessData::size
virtual int size() const
geotess::GeoTessEnumType::ordinal
int ordinal() const
Definition GeoTessEnumType.h:148
geotess::GeoTessEnumType::name
string name() const
Definition GeoTessEnumType.h:143
geotess::GeoTessException
An exception class for all GeoTess objects.
Definition GeoTessException.h:66
geotess::GeoTessInterpolatorType
Enumeration of the interpolation algorithms supported by GeoTess including LINEAR,...
Definition GeoTessInterpolatorType.h:72
geotess::GeoTessMetaData
Basic metadata information about a GeoTessModel.
Definition GeoTessMetaData.h:97
geotess::GeoTessProfile
Abstract class that manages the radii and data values that span a single layer associated with a sing...
Definition GeoTessProfile.h:97
geotess::GeoTessProfile::getData
virtual GeoTessData ** getData()
geotess::GeoTessProfile::getNRadii
virtual int getNRadii() const
geotess::GeoTessProfile::findClosestRadiusIndex
virtual int findClosestRadiusIndex(double radius) const
geotess::GeoTessProfile::getRadius
virtual float getRadius(int i) const
geotess::GeoTessProfile::getRadiusIndex
int getRadiusIndex(double radius) const
Definition GeoTessProfile.h:472
geotess::GeoTessProfileNPoint::getDataBottom
virtual GeoTessData * getDataBottom()
Definition GeoTessProfileNPoint.h:393
geotess::GeoTessProfileNPoint::GeoTessProfileNPoint
GeoTessProfileNPoint(const vector< float > &r, vector< GeoTessData * > &d)
Definition GeoTessProfileNPoint.h:185
geotess::GeoTessProfileNPoint::GeoTessProfileNPoint
GeoTessProfileNPoint(float *r, GeoTessData **dat, int size)
Definition GeoTessProfileNPoint.h:150
geotess::GeoTessProfileNPoint::getRadiusIndex
virtual int getRadiusIndex(double radius, int jlo) const
geotess::GeoTessProfileNPoint::setRadius
virtual void setRadius(int index, float radius)
Definition GeoTessProfileNPoint.h:356
geotess::GeoTessProfileNPoint::setData
virtual void setData(int index, GeoTessData *inData)
Definition GeoTessProfileNPoint.h:362
geotess::GeoTessProfileNPoint::isNaN
bool isNaN(int nodeIndex, int attributeIndex)
Definition GeoTessProfileNPoint.h:292
geotess::GeoTessProfileNPoint::setRadii
virtual void setRadii(const vector< float > &newRadii)
Definition GeoTessProfileNPoint.h:353
geotess::GeoTessProfileNPoint::getRadiusBottom
virtual float getRadiusBottom() const
Definition GeoTessProfileNPoint.h:383
geotess::GeoTessProfileNPoint::getRadii
virtual float * getRadii()
Definition GeoTessProfileNPoint.h:315
geotess::GeoTessProfileNPoint::getData
virtual const GeoTessData & getData(int i) const
Definition GeoTessProfileNPoint.h:343
geotess::GeoTessProfileNPoint::getRadius
virtual float getRadius(int i) const
Definition GeoTessProfileNPoint.h:299
geotess::GeoTessProfileNPoint::getDataTop
virtual const GeoTessData & getDataTop() const
Definition GeoTessProfileNPoint.h:373
geotess::GeoTessProfileNPoint::getValue
virtual double getValue(int attributeIndex, int radiusIndex) const
Definition GeoTessProfileNPoint.h:400
geotess::GeoTessProfileNPoint::getInterpolationCoefficient
virtual double getInterpolationCoefficient(int index, double radius, bool allowOutOfRange) const
geotess::GeoTessProfileNPoint::getDataTop
virtual GeoTessData * getDataTop()
Definition GeoTessProfileNPoint.h:378
geotess::GeoTessProfileNPoint::getNRadii
virtual int getNRadii() const
Definition GeoTessProfileNPoint.h:304
geotess::GeoTessProfileNPoint::getData
virtual GeoTessData ** getData()
Definition GeoTessProfileNPoint.h:327
geotess::GeoTessProfileNPoint::getDataBottom
virtual const GeoTessData & getDataBottom() const
Definition GeoTessProfileNPoint.h:388
geotess::GeoTessProfileNPoint::GeoTessProfileNPoint
GeoTessProfileNPoint(float *rad, const vector< GeoTessData * > &dat)
geotess::GeoTessProfileNPoint::getNData
virtual int getNData() const
Definition GeoTessProfileNPoint.h:309
geotess::GeoTessProfileNPoint::getValueTop
virtual double getValueTop(int attributeIndex) const
Definition GeoTessProfileNPoint.h:409
geotess::GeoTessProfileNPoint::class_size
virtual int class_size() const
Definition GeoTessProfileNPoint.h:235
geotess::GeoTessProfileNPoint::class_name
static string class_name()
Definition GeoTessProfileNPoint.h:230
geotess::GeoTessProfileNPoint::getData
virtual GeoTessData * getData(int i)
Definition GeoTessProfileNPoint.h:338
geotess::GeoTessProfileNPoint::getMemory
virtual LONG_INT getMemory()
Definition GeoTessProfileNPoint.h:238
geotess::GeoTessProfileNPoint::getInterpolationCoefficient
virtual double getInterpolationCoefficient(int index, double radius) const
Definition GeoTessProfileNPoint.h:649
geotess::GeoTessProfileNPoint::getType
virtual const GeoTessProfileType & getType() const
Definition GeoTessProfileNPoint.h:265
geotess::GeoTessProfileNPoint::setData
virtual void setData(const vector< GeoTessData * > &inData)
geotess::GeoTessProfileNPoint::getRadiusTop
virtual float getRadiusTop() const
Definition GeoTessProfileNPoint.h:368
geotess::GeoTessProfileType
Enumeration of the valid Profile types, including EMPTY, THIN, CONSTANT, NPOINT and SURFACE.
Definition GeoTessProfileType.h:70
geotess::GeoTessProfileType::NPOINT
static const GeoTessProfileType NPOINT
Definition GeoTessProfileType.h:118
geotess::IFStreamAscii
Opens ascii file for read and write access.
Definition IFStreamAscii.h:81
geotess::IFStreamBinary
Opens a file for binary read and write access.
Definition IFStreamBinary.h:80
geotess
Definition ArrayReuse.h:55