USGS

Isis 3.0 Object Programmers' Reference

Home

EquatorialCylindricalShape.cpp
2 
3 #include <algorithm>
4 #include <cfloat>
5 #include <string>
6 #include <iomanip>
7 #include <cmath>
8 #include <vector>
9 
10 #include <SpiceUsr.h>
11 #include <SpiceZfc.h>
12 #include <SpiceZmc.h>
13 
14 #include "Cube.h"
15 #include "IException.h"
16 #include "Latitude.h"
17 #include "Longitude.h"
18 #include "NaifStatus.h"
19 #include "SpecialPixel.h"
20 #include "SurfacePoint.h"
21 #include "Table.h"
22 
23 #define MAX(x,y) (((x) > (y)) ? (x) : (y))
24 
25 namespace Isis {
32  DemShape(target, pvl) {
33  setName("EquatorialCylindricalShape");
34 
35  m_minRadius = NULL;
36  m_maxRadius = NULL;
37 
38  // Read in the min/max radius of the DEM file and the Scale of the DEM
39  // file in pixels/degree
40  if (!demCube()->hasTable("ShapeModelStatistics")) {
41  QString msg = "The input cube references a ShapeModel that has "
42  "not been updated for the new ray tracing algorithm. All DEM "
43  "files must now be padded at the poles and contain a "
44  "ShapeModelStatistics table defining their minimum and maximum "
45  "radii values. The demprep program should be used to prepare the "
46  "DEM before you can run this program. There is more information "
47  "available in the documentation of the demprep program.";
49  }
50 
51  // Table table("ShapeModelStatistics", demCubeFile(), *demCube()->label()));
52  Table table("ShapeModelStatistics", demCube()->fileName(), *demCube()->label());
53 
54  // Find minimum and maximum radius
55  m_minRadius = new Distance(table[0]["MinimumRadius"], Distance::Kilometers);
56  m_maxRadius = new Distance(table[0]["MaximumRadius"], Distance::Kilometers);
57  }
58 
59 
64 
65  delete m_minRadius;
66  m_minRadius = NULL;
67 
68  delete m_minRadius;
69  m_maxRadius = NULL;
70  }
71 
72 
77  std::vector<double> observerPos, std::vector<double> lookDirection) {
78 
79  // try to intersect the target body ellipsoid as a first approximation
80  // for the iterative DEM intersection method (this method is in the ShapeModel base class).
81  // If the Ellipsoid intersection fails, we're done
82  if (!DemShape::intersectSurface(observerPos, lookDirection)) {
83 
85  return hasIntersection();
86 
87  SpiceDouble a = targetRadii()[0].kilometers();
88 
89  double plen=0.0;
90  SpiceDouble plat, plon, pradius;
91  SpiceDouble pB[3];
92  double maxRadiusMetersSquared =
94 
95  double cmin = cos((90.0 - 1.0 / (2.0*demScale())) * DEG2RAD);
96 
97  // Separate iteration algorithms are used for different projections -
98  // use this iteration for equatorial cylindrical type projections that
99  // failed to find an intersection with the DemShape method.
100  int maxit = 100;
101  int it = 0;
102  bool done = false;
103 
104  // Normalize the look vector
105  SpiceDouble ulookB[3];
106  ulookB[0] = lookDirection[0];
107  ulookB[1] = lookDirection[1];
108  ulookB[2] = lookDirection[2];
109  vhat_c(ulookB,ulookB);
110 
111  // Calculate the limb viewing angle to see if the line of sight is
112  // pointing away from the planet
113  SpiceDouble observer[3];
114  observer[0] = observerPos[0];
115  observer[1] = observerPos[1];
116  observer[2] = observerPos[2];
117  SpiceDouble negobserver[3];
118  vminus_c(observer, negobserver);
119  double psi0 = vsep_c(negobserver, ulookB);
120  double cospsi0 = cos(psi0);
121 
122  // If psi0 is greater than 90 degrees, then reject data as looking
123  // away from the planet and no proper tangent point exists in the
124  // direction that the spacecraft is looking
125  if (psi0 > PI/2.0) {
126  setHasIntersection(false);
127  return hasIntersection();
128  }
129 
130  // Calculate the vector to the tangent point
131  SpiceDouble tvec[3];
132  double observerdist = vnorm_c(observer);
133  tvec[0] = observer[0] + observerdist*cospsi0*ulookB[0];
134  tvec[1] = observer[1] + observerdist*cospsi0*ulookB[1];
135  tvec[2] = observer[2] + observerdist*cospsi0*ulookB[2];
136  double tlen = vnorm_c(tvec);
137 
138  // Calculate distance along look vector to first and last test point
139  double radiusDiff = maxRadiusMetersSquared - tlen * tlen;
140 
141  // Make sure radiusDiff makes sense
142  if (radiusDiff >= 0.0) {
143  radiusDiff = sqrt(radiusDiff);
144  }
145  else {
146  setHasIntersection(false);
147  return hasIntersection();
148  }
149 
150  double d0 = observerdist * cospsi0 - radiusDiff;
151  double dm = observerdist * cospsi0 + radiusDiff;
152 
153  // Set the properties at the first test observation point
154  double d = d0;
155  SpiceDouble g1[3];
156  g1[0] = observer[0] + d0*ulookB[0];
157  g1[1] = observer[1] + d0*ulookB[1];
158  g1[2] = observer[2] + d0*ulookB[2];
159  double g1len = vnorm_c(g1);
160  SpiceDouble g1lat, g1lon, g1radius;
161  reclat_c(g1,&g1radius,&g1lon,&g1lat);
162  g1lat *= RAD2DEG;
163  g1lon *= RAD2DEG;
164 
165  if (g1lon < 0.0)
166  g1lon += 360.0;
167 
168  SpiceDouble negg1[3];
169  vminus_c(g1, negg1);
170  double psi1 = vsep_c(negg1, ulookB);
171 
172  // Set dalpha to be half the grid spacing for nyquist sampling
173  //double dalpha = (PI/180.0)/(2.0*p_demScale);
174  double dalpha = MAX(cos(g1lat*DEG2RAD),cmin) / (2.0*demScale()*DEG2RAD);
175 
176  // Previous Sensor version used local version of this method with lat and lon doubles. Steven said
177  // it didn't make a significant difference in speed.
178  double r1 = (localRadius(Latitude(g1lat, Angle::Degrees),
179  Longitude(g1lon, Angle::Degrees))).kilometers();
180 
181  if (Isis::IsSpecial(r1)) {
182  setHasIntersection(false);
183  return hasIntersection();
184  }
185 
186  // Set the tolerance to a fraction of the equatorial radius, a
187  double tolerance = 3E-8 * a;
188 
189  // Main iteration loop
190  // Loop from g1 to gm stepping by angles of dalpha until intersection is found
191  while (!done) {
192 
193  if (d > dm) {
194  setHasIntersection(false);
195  return hasIntersection();
196  }
197 
198  it = 0;
199 
200  // Calculate the angle between the look vector and the planet radius at the current
201  // test point
202  double psi2 = psi1 + dalpha;
203 
204  // Calculate the step size
205  double dd = g1len * sin(dalpha) / sin(PI-psi2);
206 
207  // JAA: If we are moving along the vector at a smaller increment than the pixel
208  // tolerance we will be in an infinite loop. The infinite loop is elimnated by
209  // this test. Now the algorithm produces a jagged limb in phocube. This may
210  // be a function of the very low resolution of the Vesta DEM and could
211  // improve in the future
212  if (dd < tolerance) {
213  setHasIntersection(false);
214  return hasIntersection();
215  }
216 
217  // Calculate the vector to the current test point from the planet center
218  d = d + dd;
219  SpiceDouble g2[3];
220  g2[0] = observer[0] + d * ulookB[0];
221  g2[1] = observer[1] + d * ulookB[1];
222  g2[2] = observer[2] + d * ulookB[2];
223  double g2len = vnorm_c(g2);
224 
225  // Determine lat,lon,radius at this point
226  SpiceDouble g2lat, g2lon, g2radius;
227  reclat_c(g2,&g2radius,&g2lon,&g2lat);
228 
229  g2lat *= RAD2DEG;
230  g2lon *= RAD2DEG;
231 
232  if (g2lon < 0.0)
233  g2lon += 360.0;
234 
235  // Previous Sensor version used local version of this method with lat and lon doubles to save
236  // According to Steven the savings was negligible.
237  double r2 = (localRadius(Latitude(g2lat, Angle::Degrees),
238  Longitude(g2lon, Angle::Degrees))).kilometers();
239 
240 
241  if (Isis::IsSpecial(r2)) {
242  setHasIntersection(false);
243  return hasIntersection();
244  }
245 
246  // Test for intersection
247  if (r2 > g2len) {
248  // An intersection has occurred. Interpolate between g1 and g2 to get the
249  // lat,lon of the intersect point.
250 
251  // If g1 and g2 straddle a hill, then we may need to iterate a few times
252  // until we are on the linear segment.
253  while (it < maxit && !done) {
254  // Calculate the fractional distance "v" to move along the look vector
255  // to the intersection point. Check to see if there was a convergence
256  // of the solution and the tolerance was too small to detect it.
257  double palt;
258  if ((g2len*r1/r2 - g1len) == 0.0) {
259  setHasIntersection(true);
260  plen = pradius;
261  palt = 0.0;
262  done = true;
263  }
264  else {
265  double v = (r1-g1len) / (g2len*r1/r2 - g1len);
266  pB[0] = g1[0] + v * dd * ulookB[0];
267  pB[1] = g1[1] + v * dd * ulookB[1];
268  pB[2] = g1[2] + v * dd * ulookB[2];
269  plen = vnorm_c(pB);
270  reclat_c(pB,&pradius,&plon,&plat);
271  plat *= RAD2DEG;
272  plon *= RAD2DEG;
273  if (plon < 0.0)
274  plon += 360.0;
275 
276  if (plon > 360.0)
277  plon -= 360.0;
278 
279  // Previous Sensor version used local version of this method with lat and lon doubles. ..Why Jeff???
280  pradius = (localRadius(Latitude(plat, Angle::Degrees),
281  Longitude(plon, Angle::Degrees))).kilometers();
282 
283  if (Isis::IsSpecial(pradius)) {
284  setHasIntersection(false);
285  return hasIntersection();
286  }
287  palt = plen - pradius;
288 
289  // The altitude relative to surface is +ve at the observation point,
290  // so reset g1=p and r1=pradius
291  if (palt > tolerance) {
292  it = it + 1;
293  g1[0] = pB[0];
294  g1[1] = pB[1];
295  g1[2] = pB[2];
296  g1len = plen;
297  r1 = pradius;
298  dd = dd * (1.0 - v);
299 
300  // The altitude relative to surface -ve at the observation point,
301  // so reset g2=p and r2=pradius
302  }
303  else if (palt < -tolerance) {
304  it = it + 1;
305  g2[0] = pB[0];
306  g2[1] = pB[1];
307  g2[2] = pB[2];
308  g2len = plen;
309  r2 = pradius;
310  dd = dd * v;
311 
312  // We are within the tolerance, so the solution has converged
313  }
314  else {
315  setHasIntersection(true);
316  plen = pradius;
317  palt = 0.0;
318  done = true;
319  }
320  }
321  if (!done && it >= maxit) {
322  setHasIntersection(false);
323  return hasIntersection();
324  }
325  }
326  }
327 
328  g1[0] = g2[0];
329  g1[1] = g2[1];
330  g1[2] = g2[2];
331  g1len = g2len;
332  r1 = r2;
333  psi1 = psi2;
334 
335  // TODO: Examine how dalpha is computed at the limb for Vesta. It
336  // appears that eventually it gets really small and causes the loop
337  // to be neverending. Of course this could be happening for other
338  // limb images and we just have never tested the code. For example,
339  // Messenger with a DEM could cause the problem. As a result JAA
340  // put in a test (above) for dd being smaller than the pixel
341  // convergence tolerance. If so the loop exits without an
342  // intersection
343  dalpha = MAX(cos(g2lat*DEG2RAD),cmin) / (2.0*demScale()*DEG2RAD);
344  } // end while
345 
346  SpiceDouble intersectionPoint[3];
347  SpiceBoolean found;
348  surfpt_c(&observerPos[0], &lookDirection[0], plen, plen, plen,
349  intersectionPoint, &found);
350 
351  surfaceIntersection()->FromNaifArray(intersectionPoint);
352 
353  if (!found) {
354  setHasIntersection(false);
355  return hasIntersection();
356  }
357  else {
358  return hasIntersection();
359  }
360 
361  }
362 
363  setHasIntersection(true);
364  return hasIntersection();
365  }
366  // Do nothing since the DEM intersection was already successful
367 }
368