Isis 3.0 Programmer Reference
Back | Home
DemShape.cpp
1 #include "DemShape.h"
2 
3 // Qt third party includes
4 #include <QDebug>
5 #include <QVector>
6 
7 // c standard library third party includes
8 #include <algorithm>
9 #include <cfloat>
10 #include <cmath>
11 #include <iomanip>
12 #include <string>
13 #include <vector>
14 
15 // naif third party includes
16 #include <SpiceUsr.h>
17 #include <SpiceZfc.h>
18 #include <SpiceZmc.h>
19 
20 #include "Cube.h"
21 #include "CubeManager.h"
22 #include "Distance.h"
23 #include "EllipsoidShape.h"
24 //#include "Geometry3D.h"
25 #include "IException.h"
26 #include "Interpolator.h"
27 #include "Latitude.h"
28 //#include "LinearAlgebra.h"
29 #include "Longitude.h"
30 #include "NaifStatus.h"
31 #include "Portal.h"
32 #include "Projection.h"
33 #include "Pvl.h"
34 #include "Spice.h"
35 #include "SurfacePoint.h"
36 #include "Table.h"
37 #include "Target.h"
39 
40 using namespace std;
41 
42 namespace Isis {
48  DemShape::DemShape() : ShapeModel() {
49  setName("DemShape");
50  m_demProj = NULL;
51  m_demCube = NULL;
52  m_interp = NULL;
53  m_portal = NULL;
54  }
55 
56 
65  DemShape::DemShape(Target *target, Pvl &pvl) : ShapeModel(target) {
66  setName("DemShape");
67  m_demProj = NULL;
68  m_demCube = NULL;
69  m_interp = NULL;
70  m_portal = NULL;
71 
72  PvlGroup &kernels = pvl.findGroup("Kernels", Pvl::Traverse);
73 
74  QString demCubeFile;
75  if (kernels.hasKeyword("ElevationModel")) {
76  demCubeFile = (QString) kernels["ElevationModel"];
77  }
78  else if(kernels.hasKeyword("ShapeModel")) {
79  demCubeFile = (QString) kernels["ShapeModel"];
80  }
81 
82  m_demCube = CubeManager::Open(demCubeFile);
83 
84  // This caching algorithm works much better for DEMs than the default,
85  // regional. This is because the uniqueIOCachingAlgorithm keeps track
86  // of a history, which for something that isn't linearly processing a
87  // cube is worth it. The regional caching algorithm tosses out results
88  // from iteration 1 of setlookdirection (first algorithm) at iteration
89  // 4 and the next setimage has to re-read the data.
92  m_interp = new Interpolator(Interpolator::BiLinearType);
96 
97  // Read in the Scale of the DEM file in pixels/degree
98  const PvlGroup &mapgrp = m_demCube->label()->findGroup("Mapping", Pvl::Traverse);
99 
100  // Save map scale in pixels per degree
101  m_pixPerDegree = (double) mapgrp["Scale"];
102  }
103 
104 
107  m_demProj = NULL;
108 
109  // We do not have ownership of p_demCube
110  m_demCube = NULL;
111 
112  delete m_interp;
113  m_interp = NULL;
114 
115  delete m_portal;
116  m_portal = NULL;
117  }
118 
119 
137  bool DemShape::intersectSurface(vector<double> observerPos,
138  vector<double> lookDirection) {
139  // try to intersect the target body ellipsoid as a first approximation
140  // for the iterative DEM intersection method
141  // (this method is in the ShapeModel base class)
142 
143  bool ellipseIntersected = intersectEllipsoid(observerPos, lookDirection);
144  if (!ellipseIntersected) {
145  return false;
146  }
147 
148  double tol = resolution()/100; // 1/100 of a pixel
149  static const int maxit = 100;
150  int it = 1;
151  double dX, dY, dZ, dist2;
152  bool done = false;
153 
154  // latitude, longitude in Decimal Degrees
155  double latDD, lonDD;
156 
157  // in each iteration, the current surface intersect point is saved for
158  // comparison with the new, updated surface intersect point
159  SpiceDouble currentIntersectPt[3];
160  SpiceDouble newIntersectPt[3];
161 
162  // initialize updated surface intersection point to the ellipsoid
163  // intersection point coordinates
164  newIntersectPt[0] = surfaceIntersection()->GetX().kilometers();
165  newIntersectPt[1] = surfaceIntersection()->GetY().kilometers();
166  newIntersectPt[2] = surfaceIntersection()->GetZ().kilometers();
167 
168  double tol2 = tol * tol;
169 
171  while (!done) {
172 
173  if (it > maxit) {
174  setHasIntersection(false);
175  done = true;
176  continue;
177  }
178 
179  // The lat/lon calculations are done here by hand for speed & efficiency
180  // With doing it in the SurfacePoint class using p_surfacePoint, there
181  // is a 24% slowdown (which is significant in this very tightly looped call).
182  double t = newIntersectPt[0] * newIntersectPt[0] +
183  newIntersectPt[1] * newIntersectPt[1];
184 
185  latDD = atan2(newIntersectPt[2], sqrt(t)) * RAD2DEG;
186  lonDD = atan2(newIntersectPt[1], newIntersectPt[0]) * RAD2DEG;
187 
188  if (lonDD < 0) {
189  lonDD += 360;
190  }
191 
192  // Previous Sensor version used local version of this method with lat and lon doubles.
193  // Steven made the change to improve speed. He said the difference was negilgible.
194  Distance radiusKm = localRadius(Latitude(latDD, Angle::Degrees),
195  Longitude(lonDD, Angle::Degrees));
196 
197  if (Isis::IsSpecial(radiusKm.kilometers())) {
198  setHasIntersection(false);
199  return false;
200  }
201 
202  // save current surface intersect point for comparison with new, updated
203  // surface intersect point
204  memcpy(currentIntersectPt, newIntersectPt, 3 * sizeof(double));
205 
206  double r = radiusKm.kilometers();
207  bool status;
208  surfpt_c((SpiceDouble *) &observerPos[0], &lookDirection[0], r, r, r, newIntersectPt,
209  (SpiceBoolean*) &status);
210 
211  // LinearAlgebra::Vector point = LinearAlgebra::vector(observerPos[0],
212  // observerPos[1],
213  // observerPos[2]);
214  // LinearAlgebra::Vector direction = LinearAlgebra::vector(lookDirection[0],
215  // lookDirection[1],
216  // lookDirection[2]);
217  // QList<double> ellipsoidRadii;
218  // ellipsoidRadii << r << r << r;
219  // LinearAlgebra::Vector newPt = Geometry3D::intersect(point, direction, ellipsoidRadii);
220 
221  setHasIntersection(status);
222  if (!status) {
223  return status;
224  }
225 
226  dX = currentIntersectPt[0] - newIntersectPt[0];
227  dY = currentIntersectPt[1] - newIntersectPt[1];
228  dZ = currentIntersectPt[2] - newIntersectPt[2];
229  dist2 = (dX*dX + dY*dY + dZ*dZ) * 1000 * 1000;
230 
231  // Now recompute tolerance at updated surface point and recheck
232  if (dist2 < tol2) {
233  surfaceIntersection()->FromNaifArray(newIntersectPt);
234  tol = resolution() / 100.0;
235  tol2 = tol * tol;
236  if (dist2 < tol2) {
237  setHasIntersection(true);
238  done = true;
239  }
240  }
241 
242  it ++;
243  } // end of while loop
245 
246  return hasIntersection();
247  }
248 
249 
259 
260  Distance distance=Distance();
261 
262  if (lat.isValid() && lon.isValid()) {
264 
265  // The next if statement attempts to do the same as the previous one, but not as well so
266  // it was replaced.
267  // if (!m_demProj->IsGood())
268  // return Distance();
269 
271 
273 
275  m_demProj->WorldY(),
278  }
279 
280  return distance;
281  }
282 
283 
290  return m_pixPerDegree;
291  }
292 
293 
300  }
301 
302 
309  return m_demCube;
310  }
311 
312 
322  bool DemShape::isDEM() const {
323  return true;
324  }
325 
326 
334 
335  std::vector<SpiceDouble> normal(3);
336  if (neighborPoints.isEmpty()) {
337  normal[0] = normal[1] = normal[2] = 0.0;
338  setNormal(normal);
339  setHasNormal(false);
340  return;
341  }
342 
343  // subtract bottom from top and left from right and store results
344  double topMinusBottom[3];
345  vsub_c(neighborPoints[0], neighborPoints[1], topMinusBottom);
346  double rightMinusLeft[3];
347  vsub_c(neighborPoints[3], neighborPoints [2], rightMinusLeft);
348 
349  // take cross product of subtraction results to get normal
350  ucrss_c(topMinusBottom, rightMinusLeft, (SpiceDouble *) &normal[0]);
351 
352  // unitize normal (and do sanity check for magnitude)
353  double mag;
354  unorm_c((SpiceDouble *) &normal[0], (SpiceDouble *) &normal[0], &mag);
355 
356  if (mag == 0.0) {
357  normal[0] = normal[1] = normal[2] = 0.0;
358  setNormal(normal);
359  setHasNormal(false);
360  return;
361  }
362  else {
363  setHasNormal(true);
364  }
365 
366  // Check to make sure that the normal is pointing outward from the planet
367  // surface. This is done by taking the dot product of the normal and
368  // any one of the unitized xyz vectors. If the normal is pointing inward,
369  // then negate it.
370  double centerLookVect[3];
371  SpiceDouble pB[3];
373  unorm_c(pB, centerLookVect, &mag);
374  double dotprod = vdot_c((SpiceDouble *) &normal[0], centerLookVect);
375  if (dotprod < 0.0) {
376  vminus_c((SpiceDouble *) &normal[0], (SpiceDouble *) &normal[0]);
377  }
378 
379  setNormal(normal);
380  }
381 
382 
389  }
390 
391 
392 }
void calculateEllipsoidalSurfaceNormal()
Calculates the ellipsoidal surface normal.
Definition: ShapeModel.cpp:81
double WorldX() const
This returns the world X coordinate provided SetGround, SetCoordinate, SetUniversalGround, or SetWorld returned with success.
Definition: Projection.cpp:510
std::vector< double > normal()
Returns the local surface normal at the current intersection point.
Definition: ShapeModel.cpp:342
double HotSample()
Returns the sample coordinate of the center pixel in the buffer for the interpolator.
Portal * m_portal
Buffer used to read from the model.
Definition: DemShape.h:100
~DemShape()
Destroys the DemShape.
Definition: DemShape.cpp:106
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition: PvlObject.h:141
void calculateSurfaceNormal()
This method calculates the surface normal of the current intersection point.
Definition: DemShape.cpp:387
PixelType pixelType() const
Definition: Cube.cpp:1355
double degrees() const
Get the angle in units of Degrees.
Definition: Angle.h:245
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.
void setNormal(const std::vector< double >)
Sets the normal for the currect intersection point.
Definition: ShapeModel.cpp:398
Buffer for containing a two dimensional section of an image.
Definition: Portal.h:52
bool hasIntersection()
Returns intersection status.
Definition: ShapeModel.cpp:309
Pvl * label() const
Returns a pointer to the IsisLabel object associated with the cube.
Definition: Cube.cpp:1298
Cube * m_demCube
The cube containing the model.
Definition: DemShape.h:97
This class is designed to encapsulate the concept of a Latitude.
Definition: Latitude.h:59
double * DoubleBuffer() const
Returns the value of the shape buffer.
Definition: Buffer.h:153
Search child objects.
Definition: PvlObject.h:170
Projection * projection()
Definition: Cube.cpp:1391
void read(Blob &blob) const
This method will read data from the specified Blob object.
Definition: Cube.cpp:686
void ToNaifArray(double naifOutput[3]) const
A naif array is a c-style array of size 3.
Unless noted otherwise, the portions of Isis written by the USGS are public domain.
Distance measurement, usually in meters.
Definition: Distance.h:47
This algorithm is designed for applications that jump around between a couple of spots in the cube wi...
Projection * m_demProj
The projection of the model.
Definition: DemShape.h:98
void setHasIntersection(bool b)
Sets the flag to indicate whether this ShapeModel has an intersection.
Definition: ShapeModel.cpp:465
This class is designed to encapsulate the concept of a Longitude.
Definition: Longitude.h:52
double HotLine()
Returns the line coordinate of the center pixel in the buffer for the interpolator.
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:208
Degrees are generally considered more human readable, 0-360 is one circle, however most math does not...
Definition: Angle.h:69
int Samples()
Returns the number of samples needed by the interpolator.
double resolution()
Convenience method to get pixel resolution (m/pix) at current intersection point. ...
Definition: ShapeModel.cpp:503
Contains multiple PvlContainers.
Definition: PvlGroup.h:57
double kilometers() const
Get the displacement in kilometers.
void calculateLocalNormal(QVector< double * > cornerNeighborPoints)
This method calculates the local surface normal of the current intersection point.
Definition: DemShape.cpp:333
bool IsSpecial(const double d)
Returns if the input pixel is special.
Definition: SpecialPixel.h:199
void addCachingAlgorithm(CubeCachingAlgorithm *)
This will add the given caching algorithm to the list of attempted caching algorithms.
Definition: Cube.cpp:1514
Container for cube-like labels.
Definition: Pvl.h:135
static Cube * Open(const QString &cubeFileName)
This method calls the method OpenCube() on the static instance.
Definition: CubeManager.h:78
void SetPosition(const double sample, const double line, const int band)
Sets the line and sample position of the buffer.
Definition: Portal.h:109
bool isValid() const
This indicates whether we have a legitimate angle stored or are in an unset, or invalid, state.
Definition: Angle.cpp:110
This class is used to create and store valid Isis3 targets.
Definition: Target.h:63
Define shapes and provide utilities for Isis3 targets.
Definition: ShapeModel.h:68
Unless noted otherwise, the portions of Isis written by the USGS are public domain.
Cube * demCube()
Returns the cube defining the shape model.
Definition: DemShape.cpp:308
int Lines()
Returns the number of lines needed by the interpolator.
Pixel interpolator.
Definition: Interpolator.h:51
double m_pixPerDegree
Scale of DEM file in pixels per degree.
Definition: DemShape.h:99
Distance localRadius(const Latitude &lat, const Longitude &lon)
Gets the radius from the DEM, if we have one.
Definition: DemShape.cpp:258
void FromNaifArray(const double naifValues[3])
A naif array is a c-style array of size 3.
double demScale()
Return the scale of the DEM shape, in pixels per degree.
Definition: DemShape.cpp:289
static void CheckErrors(bool resetNaif=true)
This method looks for any naif errors that might have occurred.
Definition: NaifStatus.cpp:43
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:432
bool intersectSurface(std::vector< double > observerPos, std::vector< double > lookDirection)
Find the intersection point with the DEM.
Definition: DemShape.cpp:137
Interpolator * m_interp
Use bilinear interpolation from dem.
Definition: DemShape.h:101
virtual void calculateDefaultNormal()
This method calculates the default normal (Ellipsoid for backwards compatability) for the DemShape...
Definition: DemShape.cpp:298
void setName(QString name)
Sets the shape name.
Definition: ShapeModel.cpp:443
The distance is being specified in meters.
Definition: Distance.h:56
void setHasNormal(bool status)
Sets the flag to indicate whether this ShapeModel has a surface normal.
Definition: ShapeModel.cpp:492
double kilometers() const
Get the distance in kilometers.
Definition: Distance.cpp:118
const double RAD2DEG(57.29577951308232087679815481)
Multiplier for converting from radians to degrees.
SurfacePoint * surfaceIntersection() const
Returns the surface intersection for this ShapeModel.
Definition: ShapeModel.cpp:299
bool hasKeyword(const QString &name) const
Check to see if a keyword exists.
DemShape()
Construct a DemShape object.
Definition: DemShape.cpp:48
Unless noted otherwise, the portions of Isis written by the USGS are public domain.
double WorldY() const
This returns the world Y coordinate provided SetGround, SetCoordinate, SetUniversalGround, or SetWorld returned with success.
Definition: Projection.cpp:530
IO Handler for Isis Cubes.
Definition: Cube.h:158
bool isDEM() const
Indicates that this shape model is from a DEM.
Definition: DemShape.cpp:322

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 ISIS Support Center
File Modified: 07/12/2023 23:17:32