Loading [MathJax]/jax/output/NativeMML/config.js
Isis 3 Programmer Reference
DemShape.cpp
1 
6 /* SPDX-License-Identifier: CC0-1.0 */
7 #include "DemShape.h"
8 
9 // Qt third party includes
10 #include <QDebug>
11 #include <QVector>
12 
13 // c standard library third party includes
14 #include <algorithm>
15 #include <cfloat>
16 #include <cmath>
17 #include <iomanip>
18 #include <string>
19 #include <vector>
20 
21 // naif third party includes
22 #include <SpiceUsr.h>
23 #include <SpiceZfc.h>
24 #include <SpiceZmc.h>
25 
26 #include "Cube.h"
27 #include "CubeManager.h"
28 #include "Distance.h"
29 #include "EllipsoidShape.h"
30 //#include "Geometry3D.h"
31 #include "IException.h"
32 #include "Interpolator.h"
33 #include "Latitude.h"
34 //#include "LinearAlgebra.h"
35 #include "Longitude.h"
36 #include "NaifStatus.h"
37 #include "Portal.h"
38 #include "Projection.h"
39 #include "Pvl.h"
40 #include "Spice.h"
41 #include "SurfacePoint.h"
42 #include "Table.h"
43 #include "Target.h"
44 #include "UniqueIOCachingAlgorithm.h"
45 
46 using namespace std;
47 
48 namespace Isis {
54  DemShape::DemShape() : ShapeModel() {
55  setName("DemShape");
56  m_demProj = NULL;
57  m_demCube = NULL;
58  m_interp = NULL;
59  m_portal = NULL;
60  }
61 
62 
71  DemShape::DemShape(Target *target, Pvl &pvl) : ShapeModel(target) {
72  setName("DemShape");
73  m_demProj = NULL;
74  m_demCube = NULL;
75  m_interp = NULL;
76  m_portal = NULL;
77 
78  PvlGroup &kernels = pvl.findGroup("Kernels", Pvl::Traverse);
79 
80  QString demCubeFile;
81  if (kernels.hasKeyword("ElevationModel")) {
82  demCubeFile = (QString) kernels["ElevationModel"];
83  }
84  else if(kernels.hasKeyword("ShapeModel")) {
85  demCubeFile = (QString) kernels["ShapeModel"];
86  }
87 
88  m_demCube = CubeManager::Open(demCubeFile);
89 
90  // This caching algorithm works much better for DEMs than the default,
91  // regional. This is because the uniqueIOCachingAlgorithm keeps track
92  // of a history, which for something that isn't linearly processing a
93  // cube is worth it. The regional caching algorithm tosses out results
94  // from iteration 1 of setlookdirection (first algorithm) at iteration
95  // 4 and the next setimage has to re-read the data.
98  m_interp = new Interpolator(Interpolator::BiLinearType);
100  m_demCube->pixelType(),
102 
103  // Read in the Scale of the DEM file in pixels/degree
104  const PvlGroup &mapgrp = m_demCube->label()->findGroup("Mapping", Pvl::Traverse);
105 
106  // Save map scale in pixels per degree
107  m_pixPerDegree = (double) mapgrp["Scale"];
108  }
109 
110 
113  m_demProj = NULL;
114 
115  // We do not have ownership of p_demCube
116  m_demCube = NULL;
117 
118  delete m_interp;
119  m_interp = NULL;
120 
121  delete m_portal;
122  m_portal = NULL;
123  }
124 
125 
143  bool DemShape::intersectSurface(vector<double> observerPos,
144  vector<double> lookDirection) {
145  // try to intersect the target body ellipsoid as a first approximation
146  // for the iterative DEM intersection method
147  // (this method is in the ShapeModel base class)
148 
149  bool ellipseIntersected = intersectEllipsoid(observerPos, lookDirection);
150  if (!ellipseIntersected) {
151  return false;
152  }
153 
154  double tol = resolution()/100; // 1/100 of a pixel
155  static const int maxit = 100;
156  int it = 1;
157  double dX, dY, dZ, dist2;
158  bool done = false;
159 
160  // latitude, longitude in Decimal Degrees
161  double latDD, lonDD;
162 
163  // in each iteration, the current surface intersect point is saved for
164  // comparison with the new, updated surface intersect point
165  SpiceDouble currentIntersectPt[3];
166  SpiceDouble newIntersectPt[3];
167 
168  // initialize updated surface intersection point to the ellipsoid
169  // intersection point coordinates
170  newIntersectPt[0] = surfaceIntersection()->GetX().kilometers();
171  newIntersectPt[1] = surfaceIntersection()->GetY().kilometers();
172  newIntersectPt[2] = surfaceIntersection()->GetZ().kilometers();
173 
174  double tol2 = tol * tol;
175 
177  while (!done) {
178 
179  if (it > maxit) {
180  setHasIntersection(false);
181  done = true;
182  continue;
183  }
184 
185  // The lat/lon calculations are done here by hand for speed & efficiency
186  // With doing it in the SurfacePoint class using p_surfacePoint, there
187  // is a 24% slowdown (which is significant in this very tightly looped call).
188  double t = newIntersectPt[0] * newIntersectPt[0] +
189  newIntersectPt[1] * newIntersectPt[1];
190 
191  latDD = atan2(newIntersectPt[2], sqrt(t)) * RAD2DEG;
192  lonDD = atan2(newIntersectPt[1], newIntersectPt[0]) * RAD2DEG;
193 
194  if (lonDD < 0) {
195  lonDD += 360;
196  }
197 
198  // Previous Sensor version used local version of this method with lat and lon doubles.
199  // Steven made the change to improve speed. He said the difference was negilgible.
200  Distance radiusKm = localRadius(Latitude(latDD, Angle::Degrees),
201  Longitude(lonDD, Angle::Degrees));
202 
203  if (Isis::IsSpecial(radiusKm.kilometers())) {
204  setHasIntersection(false);
205  return false;
206  }
207 
208  // save current surface intersect point for comparison with new, updated
209  // surface intersect point
210  memcpy(currentIntersectPt, newIntersectPt, 3 * sizeof(double));
211 
212  double r = radiusKm.kilometers();
213  bool status;
214  surfpt_c((SpiceDouble *) &observerPos[0], &lookDirection[0], r, r, r, newIntersectPt,
215  (SpiceBoolean*) &status);
216 
217  // LinearAlgebra::Vector point = LinearAlgebra::vector(observerPos[0],
218  // observerPos[1],
219  // observerPos[2]);
220  // LinearAlgebra::Vector direction = LinearAlgebra::vector(lookDirection[0],
221  // lookDirection[1],
222  // lookDirection[2]);
223  // QList<double> ellipsoidRadii;
224  // ellipsoidRadii << r << r << r;
225  // LinearAlgebra::Vector newPt = Geometry3D::intersect(point, direction, ellipsoidRadii);
226 
227  setHasIntersection(status);
228  if (!status) {
229  return status;
230  }
231 
232  dX = currentIntersectPt[0] - newIntersectPt[0];
233  dY = currentIntersectPt[1] - newIntersectPt[1];
234  dZ = currentIntersectPt[2] - newIntersectPt[2];
235  dist2 = (dX*dX + dY*dY + dZ*dZ) * 1000 * 1000;
236 
237  // Now recompute tolerance at updated surface point and recheck
238  if (dist2 < tol2) {
239  surfaceIntersection()->FromNaifArray(newIntersectPt);
240  tol = resolution() / 100.0;
241  tol2 = tol * tol;
242  if (dist2 < tol2) {
243  setHasIntersection(true);
244  done = true;
245  }
246  }
247 
248  it ++;
249  } // end of while loop
251 
252  return hasIntersection();
253  }
254 
255 
265 
266  Distance distance=Distance();
267 
268  if (lat.isValid() && lon.isValid()) {
270 
271  // The next if statement attempts to do the same as the previous one, but not as well so
272  // it was replaced.
273  // if (!m_demProj->IsGood())
274  // return Distance();
275 
277 
279 
281  m_demProj->WorldY(),
284  }
285 
286  return distance;
287  }
288 
289 
296  return m_pixPerDegree;
297  }
298 
299 
306 
307  if (!surfaceIntersection()->Valid() || !hasIntersection() ) {
308  IString msg = "A valid intersection must be defined before computing the surface normal";
309  throw IException(IException::Programmer, msg, _FILEINFO_);
310  }
311 
312  // Get the coordinates of the current surface point
313  SpiceDouble pB[3];
314  pB[0] = surfaceIntersection()->GetX().kilometers();
315  pB[1] = surfaceIntersection()->GetY().kilometers();
316  pB[2] = surfaceIntersection()->GetZ().kilometers();
317 
318  // Get the radii of the ellipsoid
319  vector<Distance> radii = targetRadii();
320  double a = radii[0].kilometers();
321  double b = radii[1].kilometers();
322  double c = radii[2].kilometers();
323 
324  vector<double> normal(3,0.);
325 
327  surfnm_c(a, b, c, pB, (SpiceDouble *) &normal[0]);
329 
330  setNormal(normal);
331  setHasNormal(true);
332 
333  }
334 
335 
336 
343  return m_demCube;
344  }
345 
346 
356  bool DemShape::isDEM() const {
357  return true;
358  }
359 
360 
368 
369  std::vector<SpiceDouble> normal(3);
370  if (neighborPoints.isEmpty()) {
371  normal[0] = normal[1] = normal[2] = 0.0;
372  setNormal(normal);
373  setHasNormal(false);
374  return;
375  }
376 
377  // subtract bottom from top and left from right and store results
378  double topMinusBottom[3];
379  vsub_c(neighborPoints[0], neighborPoints[1], topMinusBottom);
380  double rightMinusLeft[3];
381  vsub_c(neighborPoints[3], neighborPoints [2], rightMinusLeft);
382 
383  // take cross product of subtraction results to get normal
384  ucrss_c(topMinusBottom, rightMinusLeft, (SpiceDouble *) &normal[0]);
385 
386  // unitize normal (and do sanity check for magnitude)
387  double mag;
388  unorm_c((SpiceDouble *) &normal[0], (SpiceDouble *) &normal[0], &mag);
389 
390  if (mag == 0.0) {
391  normal[0] = normal[1] = normal[2] = 0.0;
392  setNormal(normal);
393  setHasNormal(false);
394  return;
395  }
396  else {
397  setHasNormal(true);
398  }
399 
400  // Check to make sure that the normal is pointing outward from the planet
401  // surface. This is done by taking the dot product of the normal and
402  // any one of the unitized xyz vectors. If the normal is pointing inward,
403  // then negate it.
404  double centerLookVect[3];
405  SpiceDouble pB[3];
407  unorm_c(pB, centerLookVect, &mag);
408  double dotprod = vdot_c((SpiceDouble *) &normal[0], centerLookVect);
409  if (dotprod < 0.0) {
410  vminus_c((SpiceDouble *) &normal[0], (SpiceDouble *) &normal[0]);
411  }
412 
413  setNormal(normal);
414  }
415 
416 
423  }
424 
425 
426 }
Isis::DemShape::DemShape
DemShape()
Construct a DemShape object.
Definition: DemShape.cpp:54
Isis::Distance::kilometers
double kilometers() const
Get the distance in kilometers.
Definition: Distance.cpp:106
Isis::Angle::Degrees
@ Degrees
Degrees are generally considered more human readable, 0-360 is one circle, however most math does not...
Definition: Angle.h:56
Isis::PvlObject::findGroup
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition: PvlObject.h:129
Isis::DemShape::demScale
double demScale()
Return the scale of the DEM shape, in pixels per degree.
Definition: DemShape.cpp:295
Isis::Portal
Buffer for containing a two dimensional section of an image.
Definition: Portal.h:36
Isis::Latitude
This class is designed to encapsulate the concept of a Latitude.
Definition: Latitude.h:51
Isis::DemShape::intersectSurface
bool intersectSurface(std::vector< double > observerPos, std::vector< double > lookDirection)
Find the intersection point with the DEM.
Definition: DemShape.cpp:143
Isis::ShapeModel::surfaceIntersection
SurfacePoint * surfaceIntersection() const
Returns the surface intersection for this ShapeModel.
Definition: ShapeModel.cpp:358
Isis::DemShape::m_demProj
Projection * m_demProj
The projection of the model.
Definition: DemShape.h:96
Isis::Cube::read
void read(Blob &blob, const std::vector< PvlKeyword > keywords=std::vector< PvlKeyword >()) const
This method will read data from the specified Blob object.
Definition: Cube.cpp:807
Isis::ShapeModel::targetRadii
std::vector< Distance > targetRadii() const
Returns the radii of the body in km.
Definition: ShapeModel.cpp:465
Isis::ShapeModel::setHasIntersection
void setHasIntersection(bool b)
Sets the flag to indicate whether this ShapeModel has an intersection.
Definition: ShapeModel.cpp:554
Isis::Buffer::DoubleBuffer
double * DoubleBuffer() const
Returns the value of the shape buffer.
Definition: Buffer.h:138
Isis::DemShape::demCube
Cube * demCube()
Returns the cube defining the shape model.
Definition: DemShape.cpp:342
Isis::PvlContainer::hasKeyword
bool hasKeyword(const QString &name) const
Check to see if a keyword exists.
Definition: PvlContainer.cpp:159
Isis::DemShape::m_portal
Portal * m_portal
Buffer used to read from the model.
Definition: DemShape.h:98
Isis::Pvl
Container for cube-like labels.
Definition: Pvl.h:119
Isis::Interpolator::Lines
int Lines()
Returns the number of lines needed by the interpolator.
Definition: Interpolator.cpp:236
Isis::DemShape::calculateDefaultNormal
virtual void calculateDefaultNormal()
This method calculates the default normal (Ellipsoid for backwards compatability) for the DemShape.
Definition: DemShape.cpp:305
Isis::CubeManager::Open
static Cube * Open(const QString &cubeFileName)
This method calls the method OpenCube() on the static instance.
Definition: CubeManager.h:84
Isis::Projection::SetUniversalGround
virtual bool SetUniversalGround(const double coord1, const double coord2)
This method is used to set the lat/lon or radius/azimuth (i.e.
Definition: Projection.cpp:417
Isis::NaifStatus::CheckErrors
static void CheckErrors(bool resetNaif=true)
This method looks for any naif errors that might have occurred.
Definition: NaifStatus.cpp:28
Isis::ShapeModel::setNormal
void setNormal(const std::vector< double >)
Sets the normal for the currect intersection point.
Definition: ShapeModel.cpp:487
Isis::IsSpecial
bool IsSpecial(const double d)
Returns if the input pixel is special.
Definition: SpecialPixel.h:197
Isis::Interpolator::HotSample
double HotSample()
Returns the sample coordinate of the center pixel in the buffer for the interpolator.
Definition: Interpolator.cpp:265
Isis::DemShape::calculateLocalNormal
void calculateLocalNormal(QVector< double * > cornerNeighborPoints)
This method calculates the local surface normal of the current intersection point.
Definition: DemShape.cpp:367
Isis::Distance
Distance measurement, usually in meters.
Definition: Distance.h:34
Isis::Longitude
This class is designed to encapsulate the concept of a Longitude.
Definition: Longitude.h:40
Isis::PvlObject::Traverse
@ Traverse
Search child objects.
Definition: PvlObject.h:158
Isis::DemShape::localRadius
Distance localRadius(const Latitude &lat, const Longitude &lon)
Gets the radius from the DEM, if we have one.
Definition: DemShape.cpp:264
Isis::Projection::WorldY
virtual double WorldY() const
This returns the world Y coordinate provided SetGround, SetCoordinate, SetUniversalGround,...
Definition: Projection.cpp:544
Isis::ShapeModel::setName
void setName(QString name)
Sets the shape name.
Definition: ShapeModel.cpp:532
Isis::ShapeModel::resolution
double resolution()
Convenience method to get pixel resolution (m/pix) at current intersection point.
Definition: ShapeModel.cpp:592
Isis::DemShape::~DemShape
~DemShape()
Destroys the DemShape.
Definition: DemShape.cpp:112
Isis::Distance::Meters
@ Meters
The distance is being specified in meters.
Definition: Distance.h:43
Isis::DemShape::m_interp
Interpolator * m_interp
Use bilinear interpolation from dem.
Definition: DemShape.h:99
Isis::PvlGroup
Contains multiple PvlContainers.
Definition: PvlGroup.h:41
Isis::Cube::addCachingAlgorithm
void addCachingAlgorithm(CubeCachingAlgorithm *)
This will add the given caching algorithm to the list of attempted caching algorithms.
Definition: Cube.cpp:1922
Isis::DemShape::m_pixPerDegree
double m_pixPerDegree
Scale of DEM file in pixels per degree.
Definition: DemShape.h:97
Isis::ShapeModel::intersectEllipsoid
bool intersectEllipsoid(const std::vector< double > observerPosRelativeToTarget, const std::vector< double > &observerLookVectorToTarget)
Finds the intersection point on the ellipsoid model using the given position of the observer (spacecr...
Definition: ShapeModel.cpp:267
Isis::Projection::WorldX
virtual double WorldX() const
This returns the world X coordinate provided SetGround, SetCoordinate, SetUniversalGround,...
Definition: Projection.cpp:524
Isis::SurfacePoint::FromNaifArray
void FromNaifArray(const double naifValues[3])
A naif array is a c-style array of size 3.
Definition: SurfacePoint.cpp:891
Isis::Cube
IO Handler for Isis Cubes.
Definition: Cube.h:167
Isis::Interpolator
Pixel interpolator.
Definition: Interpolator.h:34
Isis::IException
Isis exception class.
Definition: IException.h:91
Isis::ShapeModel::setHasNormal
void setHasNormal(bool status)
Sets the flag to indicate whether this ShapeModel has a surface normal.
Definition: ShapeModel.cpp:581
Isis::Displacement::kilometers
double kilometers() const
Get the displacement in kilometers.
Definition: Displacement.cpp:94
Isis::IException::Programmer
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition: IException.h:146
std
Namespace for the standard library.
Isis::ShapeModel
Define shapes and provide utilities for Isis targets.
Definition: ShapeModel.h:62
Isis::Angle::isValid
bool isValid() const
This indicates whether we have a legitimate angle stored or are in an unset, or invalid,...
Definition: Angle.cpp:95
Isis::Cube::pixelType
PixelType pixelType() const
Definition: Cube.cpp:1758
Isis::Interpolator::Interpolate
double Interpolate(const double isamp, const double iline, const double buf[])
Performs an interpolation on the data according to the parameters set in the constructor.
Definition: Interpolator.cpp:56
Isis::Cube::label
Pvl * label() const
Returns a pointer to the IsisLabel object associated with the cube.
Definition: Cube.cpp:1701
Isis::SurfacePoint::ToNaifArray
void ToNaifArray(double naifOutput[3]) const
A naif array is a c-style array of size 3.
Definition: SurfacePoint.cpp:870
Isis::UniqueIOCachingAlgorithm
This algorithm is designed for applications that jump around between a couple of spots in the cube wi...
Definition: UniqueIOCachingAlgorithm.h:30
Isis::Angle::degrees
double degrees() const
Get the angle in units of Degrees.
Definition: Angle.h:232
Isis::ShapeModel::hasIntersection
bool hasIntersection()
Returns intersection status.
Definition: ShapeModel.cpp:368
Isis::Interpolator::Samples
int Samples()
Returns the number of samples needed by the interpolator.
Definition: Interpolator.cpp:208
Isis::IString
Adds specific functionality to C++ strings.
Definition: IString.h:165
QVector
This is free and unencumbered software released into the public domain.
Definition: Calculator.h:18
Isis::Interpolator::HotLine
double HotLine()
Returns the line coordinate of the center pixel in the buffer for the interpolator.
Definition: Interpolator.cpp:297
Isis::Portal::SetPosition
void SetPosition(const double sample, const double line, const int band)
Sets the line and sample position of the buffer.
Definition: Portal.h:93
Isis::Target
This class is used to create and store valid Isis targets.
Definition: Target.h:63
Isis::RAD2DEG
const double RAD2DEG
Multiplier for converting from radians to degrees.
Definition: Constants.h:44
Isis::Cube::projection
Projection * projection()
Definition: Cube.cpp:1794
Isis::DemShape::isDEM
bool isDEM() const
Indicates that this shape model is from a DEM.
Definition: DemShape.cpp:356
Isis::DemShape::calculateSurfaceNormal
void calculateSurfaceNormal()
This method calculates the surface normal of the current intersection point.
Definition: DemShape.cpp:421
Isis
This is free and unencumbered software released into the public domain.
Definition: Apollo.h:16
Isis::DemShape::m_demCube
Cube * m_demCube
The cube containing the model.
Definition: DemShape.h:95
Isis::ShapeModel::normal
std::vector< double > normal()
Returns the local surface normal at the current intersection point.
Definition: ShapeModel.cpp:401

U.S. Department of the Interior | U.S. Geological Survey
ISIS | Privacy & Disclaimers | Astrogeology Research Program
To contact us, please post comments and questions on the USGS Astrogeology Discussion Board
To report a bug, or suggest a feature go to: ISIS Github
File Modified: 03/21/2022 06:55:49