Isis 3 Programmer Reference
ImagePolygon.cpp
Go to the documentation of this file.
1 
24 #include "IsisDebug.h"
25 
26 #include <string>
27 #include <iostream>
28 #include <vector>
29 
30 #include <QDebug>
31 
32 #include <geos/geom/Geometry.h>
33 #include <geos/geom/Polygon.h>
34 #include <geos/geom/CoordinateArraySequence.h>
35 #include <geos/algorithm/LineIntersector.h>
36 #include <geos/util/IllegalArgumentException.h>
37 #include <geos/util/TopologyException.h>
38 #include <geos/util/GEOSException.h>
39 #include <geos/io/WKTReader.h>
40 #include <geos/io/WKTWriter.h>
41 #include <geos/operation/distance/DistanceOp.h>
42 
43 #include "ImagePolygon.h"
44 #include "IString.h"
45 #include "SpecialPixel.h"
46 #include "PolygonTools.h"
47 
48 using namespace std;
49 
50 namespace Isis {
51 
56  ImagePolygon::ImagePolygon() : Blob("Footprint", "Polygon") {
57  p_polygons = NULL;
58 
59  p_cube = NULL;
60 
61  m_leftCoord = NULL;
62  m_rightCoord = NULL;
63  m_topCoord = NULL;
64  m_botCoord = NULL;
65 
66  p_cubeStartSamp = 1;
67  p_cubeStartLine = 1;
68 
69  p_emission = 180.0;
70  p_incidence = 180.0;
71 
72  p_subpixelAccuracy = 50; //An accuracte and quick number
73 
74  p_ellipsoid = false;
75  }
76 
77 
80  delete p_polygons;
81  p_polygons = NULL;
82 
83  p_cube = NULL;
84 
85  delete m_leftCoord;
86  m_leftCoord = NULL;
87 
88  delete m_rightCoord;
89  m_rightCoord = NULL;
90 
91  delete m_topCoord;
92  m_topCoord = NULL;
93 
94  delete m_botCoord;
95  m_botCoord = NULL;
96  }
97 
98 
114  Camera * ImagePolygon::initCube(Cube &cube, int ss, int sl,
115  int ns, int nl, int band) {
116  p_gMap = new UniversalGroundMap(cube);
117  p_gMap->SetBand(band);
118 
119  p_cube = &cube;
120 
121  Camera *cam = NULL;
122  p_isProjected = false;
123 
124  try {
125  cam = cube.camera();
126  }
127  catch(IException &camError) {
128  try {
129  cube.projection();
130  p_isProjected = true;
131  }
132  catch(IException &projError) {
133  QString msg = "Can not create polygon, ";
134  msg += "cube [" + cube.fileName();
135  msg += "] is not a camera or map projection";
136 
137  IException polyError(IException::User, msg, _FILEINFO_);
138 
139  polyError.append(camError);
140  polyError.append(projError);
141 
142  throw polyError;
143  }
144  }
145  if (cam != NULL) p_isProjected = cam->HasProjection();
146 
147  // Create brick for use in SetImage
148  p_brick = new Brick(1, 1, 1, cube.pixelType());
149 
150  //------------------------------------------------------------------------
151  // Save cube number of samples and lines for later use.
152  //------------------------------------------------------------------------
153  p_cubeSamps = cube.sampleCount();
154  p_cubeLines = cube.lineCount();
155 
156  if (ns != 0) {
157  p_cubeSamps = std::min(p_cubeSamps, ss + ns);
158  }
159 
160  if (nl != 0) {
161  p_cubeLines = std::min(p_cubeLines, sl + nl);
162  }
163 
164  p_cubeStartSamp = ss;
165  p_cubeStartLine = sl;
166 
167  if (p_ellipsoid && IsLimb() && p_gMap->Camera()) {
168  try {
170  }
171  catch(IException &) {
172  std::string msg = "Cannot use an ellipsoid shape model";
173  msg += " on a limb image without a camera.";
175  }
176  }
177 
178  return cam;
179  }
180 
181 
210  void ImagePolygon::Create(Cube &cube, int sinc, int linc,
211  int ss, int sl, int ns, int nl, int band,
212  bool increasePrecision) {
213 
214  Camera *cam = NULL;
215 
216  if (sinc < 1 || linc < 1) {
217  std::string msg = "Sample and line increments must be 1 or greater";
219  }
220 
221  cam = initCube(cube, ss, sl, ns, nl, band);
222 
223  // Reduce the increment size to find a valid polygon
224  bool polygonGenerated = false;
225  while (!polygonGenerated) {
226  try {
227  p_sampinc = sinc;
228  p_lineinc = linc;
229 
230  p_pts = NULL;
231  p_pts = new geos::geom::CoordinateArraySequence();
232 
233  WalkPoly();
234 
235  polygonGenerated = true;
236  }
237  catch (IException &e) {
238  delete p_pts;
239 
240  sinc = sinc * 2 / 3;
241  linc = linc * 2 / 3;
242 
243  if (increasePrecision && (sinc > 1 || linc > 1)) {
244  }
245  else {
246  e.print(); // This should be a NAIF error
247  QString msg = "Cannot find polygon for image "
248  "[" + cube.fileName() + "]: ";
249  msg += increasePrecision ? "Cannot increase precision any further" :
250  "The increment/step size might be too large";
252  }
253  }
254  }
255 
256  /*------------------------------------------------------------------------
257  / If image contains 0/360 boundary, the polygon needs to be split up
258  / into multi polygons.
259  /-----------------------------------------------------------------------*/
260  if (cam) {
261  Pvl defaultMap;
262  cam->BasicMapping(defaultMap);
263  }
264 
265  // Create the polygon, fixing if needed
266  Fix360Poly();
267 
268  if (p_brick != 0) delete p_brick;
269 
270  if (p_gMap->Camera())
272  }
273 
274 
291  geos::geom::Coordinate ImagePolygon::FindNextPoint(geos::geom::Coordinate *currentPoint,
292  geos::geom::Coordinate lastPoint,
293  int recursionDepth) {
294  double x = lastPoint.x - currentPoint->x;
295  double y = lastPoint.y - currentPoint->y;
296  geos::geom::Coordinate result;
297 
298  // Check to see if depth is too deep (walked all the way around and found nothing)
299  if (recursionDepth > 6) {
300  return *currentPoint;
301  }
302 
303  // Check and walk in appropriate direction
304  if (x == 0.0 && y == 0.0) { // Find the starting point on an image
305  for (int line = -1 * p_lineinc; line <= 1 * p_lineinc; line += p_lineinc) {
306  for (int samp = -1 * p_sampinc; samp <= 1 * p_sampinc; samp += p_sampinc) {
307  double s = currentPoint->x + samp;
308  double l = currentPoint->y + line;
309  // Try the next left hand rule point if (s,l) does not produce a
310  // lat/long or it is not on the image.
311  if (!InsideImage(s, l) || !SetImage(s, l)) {
312  geos::geom::Coordinate next(s, l);
313  return FindNextPoint(currentPoint, next);
314  }
315  }
316  }
317 
318  std::string msg = "Unable to create image footprint. Starting point is not on the edge of the image.";
320  }
321  else if (x < 0 && y < 0) { // current is top left
322  geos::geom::Coordinate next(currentPoint->x, currentPoint->y - 1 * p_lineinc);
323  MoveBackInsideImage(next.x, next.y, 0, -p_lineinc);
324  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
325  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
326  }
327  else {
328  result = FindBestPoint(currentPoint, next, lastPoint);
329  }
330  }
331  else if (x == 0.0 && y < 0) { // current is top
332  geos::geom::Coordinate next(currentPoint->x + 1 * p_sampinc, currentPoint->y - 1 * p_lineinc);
333  MoveBackInsideImage(next.x, next.y, p_sampinc, -p_lineinc);
334  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
335  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
336  }
337  else {
338  result = FindBestPoint(currentPoint, next, lastPoint);
339  }
340  }
341  else if (x > 0 && y < 0) { // current is top right
342  geos::geom::Coordinate next(currentPoint->x + 1 * p_sampinc, currentPoint->y);
343  MoveBackInsideImage(next.x, next.y, p_sampinc, 0);
344  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
345  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
346  }
347  else {
348  result = FindBestPoint(currentPoint, next, lastPoint);
349  }
350  }
351  else if (x > 0 && y == 0.0) { // current is right
352  geos::geom::Coordinate next(currentPoint->x + 1 * p_sampinc, currentPoint->y + 1 * p_lineinc);
353  MoveBackInsideImage(next.x, next.y, p_sampinc, p_lineinc);
354  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
355  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
356  }
357  else {
358  result = FindBestPoint(currentPoint, next, lastPoint);
359  }
360  }
361  else if (x > 0 && y > 0) { // current is bottom right
362  geos::geom::Coordinate next(currentPoint->x, currentPoint->y + 1 * p_lineinc);
363  MoveBackInsideImage(next.x, next.y, 0, p_lineinc);
364  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
365  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
366  }
367  else {
368  result = FindBestPoint(currentPoint, next, lastPoint);
369  }
370  }
371  else if (x == 0.0 && y > 0) { // current is bottom
372  geos::geom::Coordinate next(currentPoint->x - 1 * p_sampinc, currentPoint->y + 1 * p_lineinc);
373  MoveBackInsideImage(next.x, next.y, -p_sampinc, p_lineinc);
374  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
375  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
376  }
377  else {
378  result = FindBestPoint(currentPoint, next, lastPoint);
379  }
380  }
381  else if (x < 0 && y > 0) { // current is bottom left
382  geos::geom::Coordinate next(currentPoint->x - 1 * p_sampinc, currentPoint->y);
383  MoveBackInsideImage(next.x, next.y, -p_sampinc, 0);
384  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
385  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
386  }
387  else {
388  result = FindBestPoint(currentPoint, next, lastPoint);
389  }
390  }
391  else if (x < 0 && y == 0.0) { // current is left
392  geos::geom::Coordinate next(currentPoint->x - 1 * p_sampinc, currentPoint->y - 1 * p_lineinc);
393  MoveBackInsideImage(next.x, next.y, -p_sampinc, -p_lineinc);
394  if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
395  result = FindNextPoint(currentPoint, next, recursionDepth + 1);
396  }
397  else {
398  result = FindBestPoint(currentPoint, next, lastPoint);
399  }
400  }
401  else {
402  std::string msg = "Unable to create image footprint. Error walking image.";
404  }
405 
406 
407  return result;
408  }
409 
420  void ImagePolygon::MoveBackInsideImage(double &sample, double &line, double sinc, double linc) {
421  // snap to centers of pixels!
422 
423  // Starting sample to snap to
424  const double startSample = p_cubeStartSamp;
425  // Ending sample to snap to
426  const double endSample = p_cubeSamps;
427  // Starting line to snap to
428  const double startLine = p_cubeStartLine;
429  // Ending line to snap to
430  const double endLine = p_cubeLines;
431  // Original sample for this point (before sinc)
432  const double origSample = sample - sinc;
433  // Original line for this point (before linc)
434  const double origLine = line - linc;
435 
436  // We moved left off the image - snap to the edge
437  if (sample < startSample && sinc < 0) {
438  // don't snap if we started at the edge
439  if (origSample == startSample) {
440  return;
441  }
442 
443  sample = startSample;
444  }
445 
446  // We moved right off the image - snap to the edge
447  if (sample > endSample && sinc > 0) {
448  // don't snap if we started at the edge
449  if (origSample == endSample) {
450  return;
451  }
452 
453  sample = endSample;
454  }
455 
456  // We moved up off the image - snap to the edge
457  if (line < startLine && linc < 0) {
458  // don't snap if we started at the edge
459  if (origLine == startLine) {
460  return;
461  }
462 
463  line = startLine;
464  }
465 
466  // We moved down off the image - snap to the edge
467  if (line > endLine && linc > 0) {
468  // don't snap if we started at the edge
469  if (fabs(origLine - endLine) < 0.5) {
470  return;
471  }
472 
473  line = endLine;
474  }
475 
476  return;
477  }
478 
479 
488  bool ImagePolygon::InsideImage(double sample, double line) {
489  return (sample >= p_cubeStartSamp - 0.5 &&
490  line > p_cubeStartLine - 0.5 &&
491  sample <= p_cubeSamps + 0.5 &&
492  line <= p_cubeLines + 0.5);
493  }
494 
495 
502  geos::geom::Coordinate ImagePolygon::FindFirstPoint() {
503  // @todo: Brute force method, should be improved
504  for (int sample = p_cubeStartSamp; sample <= p_cubeSamps; sample++) {
505  for (int line = p_cubeStartLine; line <= p_cubeLines; line++) {
506  if (SetImage(sample, line)) {
507  // An outlier check. Make sure that the pixel we use to start
508  // constructing a polygon is not surrounded by a bunch of invalid
509  // positions.
510  geos::geom::Coordinate firstPoint(sample, line);
511  geos::geom::Coordinate lastPoint = firstPoint;
512  if (!firstPoint.equals(FindNextPoint(&firstPoint, lastPoint))) {
513  return firstPoint;
514  }
515  }
516  }
517  }
518 
519  // Check to make sure a point was found
520  std::string msg = "No lat/lon data found for image";
522  }
523 
524 
530  double result = 0.0;
531 
533  if (m_rightCoord && m_leftCoord)
534  result = m_rightCoord->x - m_leftCoord->x + 1;
535 
536  return result;
537  }
538 
539 
545  double result = 0.0;
546 
548  if (m_topCoord && m_botCoord)
549  result = m_botCoord->y - m_topCoord->y + 1;
550 
551  return result;
552  }
553 
554 
564  ASSERT(p_cube);
565 
566  for (int line = p_cubeStartLine; !m_leftCoord && line <= p_cubeLines; line++) {
567  for (int sample = p_cubeStartSamp; !m_leftCoord && sample <= p_cubeSamps; sample++) {
568  if (SetImage(sample, line)) {
569  m_leftCoord = new geos::geom::Coordinate(sample, line);
570  }
571  }
572  }
573 
574  if (m_leftCoord) {
575  for (int line = p_cubeStartLine; !m_rightCoord && line <= p_cubeLines; line++) {
576  for (int sample = p_cubeSamps; !m_rightCoord && sample >= m_leftCoord->x; sample--) {
577  if (SetImage(sample, line)) {
578  m_rightCoord = new geos::geom::Coordinate(sample, line);
579  }
580  }
581  }
582  }
583 
584  if (m_leftCoord && m_rightCoord) {
585  for (int sample = (int)m_leftCoord->x; !m_topCoord && sample <= m_rightCoord->x; sample++) {
586  for (int line = 1; !m_topCoord && line <= p_cubeLines; line++) {
587  if (SetImage(sample, line)) {
588  m_topCoord = new geos::geom::Coordinate(sample, line);
589  }
590  }
591  }
592  }
593 
594  if (m_leftCoord && m_rightCoord && m_topCoord) {
595  for (int sample = (int)m_leftCoord->x; !m_botCoord && sample <= m_rightCoord->x; sample++) {
596  for (int line = p_cube->lineCount(); !m_botCoord && line >= m_topCoord->y; line--) {
597  if (SetImage(sample, line)) {
598  m_botCoord = new geos::geom::Coordinate(sample, line);
599  }
600  }
601  }
602  }
603  }
604 
605 
614  vector<geos::geom::Coordinate> points;
615  double lat, lon, prevLat, prevLon;
616 
617  // Find the edge of the polygon
618  geos::geom::Coordinate firstPoint = FindFirstPoint();
619  points.push_back(firstPoint);
620  //************************
621  // Start walking the edge
622  //************************
623 
624  // set up for intialization
625  geos::geom::Coordinate currentPoint = firstPoint;
626  geos::geom::Coordinate lastPoint = firstPoint;
627  geos::geom::Coordinate tempPoint;
628 
629  do {
630  tempPoint = FindNextPoint(&currentPoint, lastPoint);
631  //exit(1);
632 
633 
634  // First check if the distance is within range of skipping
635  bool snapToFirstPoint = true;
636 
637  // Never needs to find firstPoint on incements of 1
638  snapToFirstPoint &= (p_sampinc != 1 && p_lineinc != 1);
639 
640  // Prevents catching the first point as the last
641  snapToFirstPoint &= (points.size() > 2);
642 
643  // This method fails for steps larger than line/sample length
644  snapToFirstPoint &= (p_sampinc < p_cubeSamps && p_lineinc < p_cubeLines);
645 
646  // Checks for appropriate distance without sqrt() call
647  int minStepSize = std::min(p_sampinc, p_lineinc);
648  snapToFirstPoint &= (DistanceSquared(&currentPoint, &firstPoint) < (minStepSize * minStepSize));
649 
650  // Searches for skipped firstPoint due to a large pixinc, by using a pixinc of 1
651  if (snapToFirstPoint) {
652  tempPoint = firstPoint;
653  }
654  // If pixinc is greater than line or sample dimention, then we could
655  // skip firstPoint. This checks for that case.
656  else if (p_sampinc > p_cubeSamps || p_lineinc > p_cubeLines) {
657  // This is not expensive because incement must be large
658  for (int pt = 0; pt < (int)points.size(); pt ++) {
659  if (points[pt].equals(tempPoint)) {
660  tempPoint = firstPoint;
661  break;
662  }
663  }
664  }
665 
666  // Failed to find the next point
667  if (tempPoint.equals(currentPoint)) {
668  geos::geom::Coordinate oldDuplicatePoint = tempPoint;
669 
670  // Init vars for first run through the loop
671  tempPoint = lastPoint;
672  lastPoint = currentPoint;
673  currentPoint = tempPoint;
674 
675  // Must be 3 (not 2) to prevent the revisit of the starting point,
676  // resulting in an infinite loop
677  if (points.size() < 3) {
678  std::string msg = "Failed to find next point in the image.";
680  }
681 
682  // remove last point from the list
683  points.pop_back();
684 
685  tempPoint = FindNextPoint(&currentPoint, lastPoint, 1);
686 
687  if (tempPoint.equals(currentPoint) || tempPoint.equals(oldDuplicatePoint)) {
688  std::string msg = "Failed to find next valid point in the image.";
690  }
691  }
692 
693 
694  // Check for triangle cycles and try to fix
695  if ((p_sampinc > 1 || p_lineinc > 1) && points.size() >= 3) {
696  if (points[points.size()-3].x == tempPoint.x &&
697  points[points.size()-3].y == tempPoint.y) {
698  // Remove the triangle from the list
699  points.pop_back();
700  points.pop_back();
701  points.pop_back();
702  // Reset the current (to be last) point
703  currentPoint = points[points.size()-1];
704  // change increment to prevent randomly bad pixels in the image
705  if (p_sampinc > 1) p_sampinc --;
706  if (p_lineinc > 1) p_lineinc --;
707  }
708 
715  if (points.size() > 250) {
716  int cycleStart = 0;
717  int cycleEnd = 0;
718 
719  for (unsigned int pt = 1; pt < points.size() && cycleStart == 0; pt ++) {
720  for (unsigned int check = pt + 1; check < points.size() && cycleStart == 0; check ++) {
721  if (points[pt] == points[check]) {
722  cycleStart = pt;
723  cycleEnd = check;
724  }
725  }
726  }
727 
728  // If a cycle was found, make it the polygon
729  if (cycleStart != 0) {
730  vector<geos::geom::Coordinate> cyclePoints;
731  for (int pt = cycleStart; pt <= cycleEnd; pt ++) {
732  cyclePoints.push_back(points[pt]);
733  }
734 
735  points = cyclePoints;
736  break;
737  }
738  }
739 
740  }
741 
742  lastPoint = currentPoint;
743  currentPoint = tempPoint;
744  points.push_back(currentPoint);
745 
746  }
747  while (!currentPoint.equals(firstPoint));
748 
749  if (points.size() <= 3) {
750  std::string msg = "Failed to find enough points on the image.";
752  }
753 
754  FindSubpixel(points);
755 
756  prevLat = 0;
757  prevLon = 0;
758  // this vector stores crossing points, where the image crosses the
759  // meridian. It stores the first coordinate of the pair in its vector
760  vector<geos::geom::Coordinate> *crossingPoints = new vector<geos::geom::Coordinate>;
761  for (unsigned int i = 0; i < points.size(); i++) {
762  geos::geom::Coordinate *temp = &(points.at(i));
763  SetImage(temp->x, temp->y);
764  lon = p_gMap->UniversalLongitude();
765  lat = p_gMap->UniversalLatitude();
766  if (abs(lon - prevLon) >= 180 && i != 0) {
767  crossingPoints->push_back(geos::geom::Coordinate(prevLon, prevLat));
768  }
769  p_pts->add(geos::geom::Coordinate(lon, lat));
770  prevLon = lon;
771  prevLat = lat;
772  }
773 
774 
775  // Checks for self-intersection and attempts to correct
776  geos::geom::CoordinateSequence *tempPts = new geos::geom::CoordinateArraySequence();
777 
778  // Gets the starting, second, second to last, and last points to check for validity
779  tempPts->add(geos::geom::Coordinate((*p_pts)[0].x, (*p_pts)[0].y));
780  tempPts->add(geos::geom::Coordinate((*p_pts)[1].x, (*p_pts)[1].y));
781  tempPts->add(geos::geom::Coordinate((*p_pts)[p_pts->size()-3].x, (*p_pts)[p_pts->size()-3].y));
782  tempPts->add(geos::geom::Coordinate((*p_pts)[p_pts->size()-2].x, (*p_pts)[p_pts->size()-2].y));
783  tempPts->add(geos::geom::Coordinate((*p_pts)[0].x, (*p_pts)[0].y));
784 
785  geos::geom::Polygon *tempPoly = globalFactory.createPolygon
786  (globalFactory.createLinearRing(tempPts), NULL);
787 
788  // Remove the last point of the sequence if it produces invalid polygons
789  if (!tempPoly->isValid()) {
790  p_pts->deleteAt(p_pts->size() - 2);
791  }
792 
793  delete tempPts;
794  tempPts = NULL;
795  // end self-intersection check
796 
797  FixPolePoly(crossingPoints);
798  delete crossingPoints;
799  crossingPoints = NULL;
800  }
801 
802 
810  void ImagePolygon::FixPolePoly(std::vector<geos::geom::Coordinate> *crossingPoints) {
811  // We currently do not support both poles in one image
812  bool hasNorthPole = false;
813  bool hasSouthPole = false;
814 
815  if (p_gMap->SetUniversalGround(90, 0)) {
816  double nPoleSample = Null;
817  double nPoleLine = Null;
818 
819  if (p_gMap->Projection()) {
820  nPoleSample = p_gMap->Projection()->WorldX();
821  nPoleLine = p_gMap->Projection()->WorldY();
822  }
823  else if (p_gMap->Camera()) {
824  nPoleSample = p_gMap->Camera()->Sample();
825  nPoleLine = p_gMap->Camera()->Line();
826  }
827 
828  if (nPoleSample >= 0.5 && nPoleLine >= 0.5 &&
829  nPoleSample <= p_cube->sampleCount() + 0.5 &&
830  nPoleLine <= p_cube->lineCount() + 0.5) {
831  hasNorthPole = true;
832  }
833  }
834 
835  if (p_gMap->SetUniversalGround(-90, 0)) {
836  double sPoleSample = Null;
837  double sPoleLine = Null;
838 
839  if (p_gMap->Projection()) {
840  sPoleSample = p_gMap->Projection()->WorldX();
841  sPoleLine = p_gMap->Projection()->WorldY();
842  }
843  else if (p_gMap->Camera()) {
844  sPoleSample = p_gMap->Camera()->Sample();
845  sPoleLine = p_gMap->Camera()->Line();
846  }
847 
848  if (sPoleSample >= 0.5 && sPoleLine >= 0.5 &&
849  sPoleSample <= p_cube->sampleCount() + 0.5 &&
850  sPoleLine <= p_cube->lineCount() + 0.5) {
851  hasSouthPole = true;
852  }
853  }
854 
855  if (hasNorthPole && hasSouthPole) {
856  std::string msg = "Unable to create image footprint because image has both poles";
858  }
859  else if (crossingPoints->size() == 0) {
860  // No crossing points
861  return;
862  }
863 
864  if (hasNorthPole) {
865  p_gMap->SetUniversalGround(90, 0);
866 
867  // If the (north) pole is settable but not within proper angles,
868  // then the polygon does not contain the (north) pole when the cube does
869  if (p_gMap->Camera() &&
871  return;
872  }
873  if (p_gMap->Camera() &&
875  return;
876  }
877  }
878  else if (hasSouthPole) {
879  p_gMap->SetUniversalGround(-90, 0);
880 
881  // If the (south) pole is settable but not within proper angles,
882  // then the polygon does not contain the (south) pole when the cube does
883  if (p_gMap->Camera() &&
885  return;
886  }
887  if (p_gMap->Camera() &&
889  return;
890  }
891  }
892 
893  geos::geom::Coordinate *closestPoint = &crossingPoints->at(0);
894  geos::geom::Coordinate *pole = NULL;
895 
896  // Setup the right pole
897  if (hasNorthPole) {
898  pole = new geos::geom::Coordinate(0, 90);
899  }
900  else if (hasSouthPole) {
901  pole = new geos::geom::Coordinate(0, -90);
902  }
903  else if (crossingPoints->size() % 2 == 1) {
904  geos::geom::Coordinate nPole(0, 90);
905  double nDist = DBL_MAX;
906  geos::geom::Coordinate sPole(0, -90);
907  double sDist = DBL_MAX;
908 
909  for (unsigned int index = 0; index < p_pts->size(); index ++) {
910  double north = DistanceSquared(&nPole, &((*p_pts)[index]));
911  if (north < nDist) {
912  nDist = north;
913  }
914  double south = DistanceSquared(&sPole, &((*p_pts)[index]));
915  if (south < sDist) {
916  sDist = south;
917  }
918  }
919 
920  if (sDist < nDist) {
921  pole = new geos::geom::Coordinate(0, -90);
922  }
923  else {
924  pole = new geos::geom::Coordinate(0, 90);
925  }
926  }
927 
928  // Skip the fix if no pole was determined
929  if (pole == NULL) {
930  return;
931  }
932 
933  // Find where the pole needs to be split
934  double closestDistance = DBL_MAX;
935  for (unsigned int i = 0; i < crossingPoints->size(); i++) {
936  geos::geom::Coordinate *temp = &crossingPoints->at(i);
937  double tempDistance;
938  // Make sure the Lat is as close to 0 as possible for a correct distance
939  if (temp->x > 180) {
940  double mod = 0.0;
941  while ((temp->x - mod) > 180) mod += 360.0;
942 
943  // Create the modified Point, create a pointer to it, and find the distance
944  geos::geom::Coordinate modPointMem = geos::geom::Coordinate(temp->x - mod, temp->y);
945  geos::geom::Coordinate *modPoint = &modPointMem;
946  tempDistance = DistanceSquared(modPoint, pole);
947  }
948  else {
949  tempDistance = DistanceSquared(temp, pole);
950  }
951  if (tempDistance < closestDistance) {
952  closestDistance = tempDistance;
953  closestPoint = temp;
954  }
955  }
956 
957  if (closestDistance == DBL_MAX) {
958  std::string msg = "Image contains a pole but did not detect a meridian crossing!";
960  }
961 
962  // split image at the pole
963  geos::geom::CoordinateSequence *new_points = new geos::geom::CoordinateArraySequence();
964  for (unsigned int i = 0; i < p_pts->size(); i++) {
965  geos::geom::Coordinate temp = p_pts->getAt(i);
966  new_points->add(temp);
967  if (temp.equals(*closestPoint)) {
968  // Setup direction
969  if (i + 1 != p_pts->size()) { // Prevent overflow exception
970  double fromLon, toLon;
971  if ((p_pts->getAt(i + 1).x - closestPoint->x) > 0) {
972  fromLon = 0.0;
973  toLon = 360.0;
974  }
975  else {
976  fromLon = 360.0;
977  toLon = 0.0;
978  }
979 
980  geos::algorithm::LineIntersector lineIntersector;
981  geos::geom::Coordinate crossingPoint;
982  geos::geom::Coordinate nPole(0.0, 90.0);
983  geos::geom::Coordinate sPole(0.0, -90.0);
984  double dist = DBL_MAX;
985 
986  for (int num = 0; num < 2 && dist > 180.0; num ++) {
987  nPole = geos::geom::Coordinate(num * 360.0, 90.0);
988  sPole = geos::geom::Coordinate(num * 360.0, -90.0);
989 
990  if (temp.x > 0.0 && p_pts->getAt(i + 1).x > 0.0) {
991  crossingPoint = geos::geom::Coordinate(p_pts->getAt(i + 1).x - 360.0 + (num * 720.0), p_pts->getAt(i + 1).y);
992  }
993  else if (temp.x < 0.0 && p_pts->getAt(i + 1).x < 0.0) { // This should never happen
994  crossingPoint = geos::geom::Coordinate(p_pts->getAt(i + 1).x + 360.0 - (num * 720.0), p_pts->getAt(i + 1).y);
995  }
996 
997  dist = std::sqrt(DistanceSquared(&temp, &crossingPoint));
998  }
999 
1000  lineIntersector.computeIntersection(nPole, sPole, temp, crossingPoint);
1001 
1002  if (lineIntersector.hasIntersection()) {
1003  const geos::geom::Coordinate &intersection = lineIntersector.getIntersection(0);
1004 
1005  // Calculate the latituded of the points along the meridian
1006  if (pole->y < intersection.y) {
1007  dist = -dist;
1008  }
1009  vector<double> lats;
1010  double maxLat = std::max(intersection.y, pole->y);
1011  double minLat = std::min(intersection.y, pole->y);
1012  for (double lat = intersection.y + dist; lat < maxLat && lat > minLat; lat += dist) {
1013  lats.push_back(lat);
1014  }
1015 
1016  // Add the new points
1017  new_points->add(geos::geom::Coordinate(fromLon, intersection.y));
1018  for (int lat = 0; lat < (int)lats.size(); lat ++) {
1019  new_points->add(geos::geom::Coordinate(fromLon, lats[lat]));
1020  }
1021  new_points->add(geos::geom::Coordinate(fromLon, pole->y));
1022  new_points->add(geos::geom::Coordinate(toLon, pole->y));
1023  for (int lat = lats.size() - 1; lat >= 0; lat --) {
1024  new_points->add(geos::geom::Coordinate(toLon, lats[lat]));
1025  }
1026  new_points->add(geos::geom::Coordinate(toLon, intersection.y));
1027  }
1028  else {
1029  std::string msg = "Image contains a pole but could not determine a meridian crossing!";
1031  }
1032 
1033  }
1034  }
1035  }
1036  delete p_pts;
1037  p_pts = new_points;
1038  delete pole;
1039  }
1040 
1041 
1054  bool ImagePolygon::SetImage(const double sample, const double line) {
1055  bool found = false;
1056  if (!p_isProjected) {
1057  found = p_gMap->SetImage(sample, line);
1058  if (!found) {
1059  return false;
1060  }
1061  else {
1062  // Check for valid emission and incidence
1063  try {
1064  if (p_gMap->Camera() &&
1066  return false;
1067  }
1068  if (p_gMap->Camera() &&
1070  return false;
1071  }
1072  }
1073  catch(IException &error) {
1074  }
1075 
1084  // Make sure we can go back to image coordinates
1085  // This is done because some camera models due to distortion, get
1086  // a lat/lon for samp/line=1:1, but entering that lat/lon does
1087  // not return samp/line =1:1. Ie. moc WA global images
1088  //double lat = p_gMap->UniversalLatitude();
1089  //double lon = p_gMap->UniversalLongitude();
1090  //return p_gMap->SetUniversalGround(lat,lon);
1091 
1092  return found;
1093  }
1094  }
1095  else {
1096  // If projected, make sure the pixel DN is valid before worrying about
1097  // geometry.
1098  p_brick->SetBasePosition((int)sample, (int)line, 1);
1099  p_cube->read(*p_brick);
1100  if (Isis::IsNullPixel((*p_brick)[0])) {
1101  return false;
1102  }
1103  else {
1104  return p_gMap->SetImage(sample, line);
1105  }
1106  }
1107  }
1108 
1109 
1117  bool convertLon = false;
1118  bool negAdjust = false;
1119  bool newCoords = false; // coordinates have been adjusted
1120  geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
1121  double lon, lat;
1122  double lonOffset = 0;
1123  double prevLon = p_pts->getAt(0).x;
1124  double prevLat = p_pts->getAt(0).y;
1125 
1126  newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
1127  double dist = 0.0;
1128  for (unsigned int i = 1; i < p_pts->getSize(); i++) {
1129  lon = p_pts->getAt(i).x;
1130  lat = p_pts->getAt(i).y;
1131  // check to see if you just crossed the Meridian
1132  if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
1133  newCoords = true;
1134  // if you were already converting then stop (crossed Meridian even number of times)
1135  if (convertLon) {
1136  convertLon = false;
1137  lonOffset = 0;
1138  }
1139  else { // Need to start converting again, deside how to adjust coordinates
1140  if ((lon - prevLon) > 0) {
1141  lonOffset = -360.;
1142  negAdjust = true;
1143  }
1144  else if ((lon - prevLon) < 0) {
1145  lonOffset = 360.;
1146  negAdjust = false;
1147  }
1148  convertLon = true;
1149  }
1150 
1151 
1152  }
1153 
1154  // Change to a minimum calculation
1155  if (newCoords && dist == 0.0) {
1156  double longitude = (lon + lonOffset) - prevLon;
1157  double latitude = lat - prevLat;
1158  dist = std::sqrt((longitude * longitude) + (latitude * latitude));
1159  }
1160 
1161  // add coord
1162  newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));
1163 
1164  // set current to old
1165  prevLon = lon;
1166  prevLat = lat;
1167  }
1168 
1169  // Nothing was done so return
1170  if (!newCoords) {
1171  geos::geom::Polygon *newPoly = globalFactory.createPolygon
1172  (globalFactory.createLinearRing(newLonLatPts), NULL);
1174  delete newLonLatPts;
1175  return;
1176  }
1177 
1178  // bisect into seperate polygons
1179  try {
1180  geos::geom::Polygon *newPoly = globalFactory.createPolygon
1181  (globalFactory.createLinearRing(newLonLatPts), NULL);
1182 
1183  geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateArraySequence();
1184  geos::geom::CoordinateSequence *pts2 = new geos::geom::CoordinateArraySequence();
1185 
1186  // Depending on direction of compensation bound accordingly
1187  //***************************************************
1188 
1189  // please verify correct if you change these values
1190  //***************************************************
1191  if (negAdjust) {
1192  pts->add(geos::geom::Coordinate(0., 90.));
1193  pts->add(geos::geom::Coordinate(-360., 90.));
1194  pts->add(geos::geom::Coordinate(-360., -90.));
1195  pts->add(geos::geom::Coordinate(0., -90.));
1196  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1197  pts->add(geos::geom::Coordinate(0.0, lat));
1198  }
1199  pts->add(geos::geom::Coordinate(0., 90.));
1200  pts2->add(geos::geom::Coordinate(0., 90.));
1201  pts2->add(geos::geom::Coordinate(360., 90.));
1202  pts2->add(geos::geom::Coordinate(360., -90.));
1203  pts2->add(geos::geom::Coordinate(0., -90.));
1204  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1205  pts2->add(geos::geom::Coordinate(0.0, lat));
1206  }
1207  pts2->add(geos::geom::Coordinate(0., 90.));
1208  }
1209  else {
1210  pts->add(geos::geom::Coordinate(360., 90.));
1211  pts->add(geos::geom::Coordinate(720., 90.));
1212  pts->add(geos::geom::Coordinate(720., -90.));
1213  pts->add(geos::geom::Coordinate(360., -90.));
1214  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1215  pts->add(geos::geom::Coordinate(360.0, lat));
1216  }
1217  pts->add(geos::geom::Coordinate(360., 90.));
1218  pts2->add(geos::geom::Coordinate(360., 90.));
1219  pts2->add(geos::geom::Coordinate(0., 90.));
1220  pts2->add(geos::geom::Coordinate(0., -90.));
1221  pts2->add(geos::geom::Coordinate(360., -90.));
1222  for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1223  pts2->add(geos::geom::Coordinate(360.0, lat));
1224  }
1225  pts2->add(geos::geom::Coordinate(360., 90.));
1226  }
1227 
1228  geos::geom::Polygon *boundaryPoly = globalFactory.createPolygon
1229  (globalFactory.createLinearRing(pts), NULL);
1230  geos::geom::Polygon *boundaryPoly2 = globalFactory.createPolygon
1231  (globalFactory.createLinearRing(pts2), NULL);
1232  /*------------------------------------------------------------------------
1233  / Intersecting the original polygon (converted coordinates) with the
1234  / boundary polygons will create the multi polygons with the converted coordinates.
1235  / These will need to be converted back to the original coordinates.
1236  /-----------------------------------------------------------------------*/
1237  geos::geom::Geometry *intersection = PolygonTools::Intersect(newPoly, boundaryPoly);
1238  geos::geom::MultiPolygon *convertPoly = PolygonTools::MakeMultiPolygon(intersection);
1239  delete intersection;
1240 
1241  intersection = PolygonTools::Intersect(newPoly, boundaryPoly2);
1242  geos::geom::MultiPolygon *convertPoly2 = PolygonTools::MakeMultiPolygon(intersection);
1243  delete intersection;
1244 
1245  /*------------------------------------------------------------------------
1246  / Adjust points created in the negative space or >360 space to be back in
1247  / the 0-360 world. This will always only need to be done on convertPoly.
1248  / Then add geometries to finalpolys.
1249  /-----------------------------------------------------------------------*/
1250  vector<geos::geom::Geometry *> *finalpolys = new vector<geos::geom::Geometry *>;
1251  geos::geom::Geometry *newGeom = NULL;
1252 
1253  for (unsigned int i = 0; i < convertPoly->getNumGeometries(); i++) {
1254  newGeom = (convertPoly->getGeometryN(i))->clone();
1255  pts = convertPoly->getGeometryN(i)->getCoordinates();
1256  geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateArraySequence();
1257  // fix the points
1258 
1259  for (unsigned int k = 0; k < pts->getSize() ; k++) {
1260  double lon = pts->getAt(k).x;
1261  double lat = pts->getAt(k).y;
1262  if (negAdjust) {
1263  lon = lon + 360;
1264  }
1265  else {
1266  lon = lon - 360;
1267  }
1268  newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);
1269 
1270  }
1271  // Add the points to polys
1272  finalpolys->push_back(globalFactory.createPolygon
1273  (globalFactory.createLinearRing(newLonLatPts), NULL));
1274  }
1275 
1276  // This loop is over polygons that will always be in 0-360 space no need to convert
1277  for (unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
1278  newGeom = (convertPoly2->getGeometryN(i))->clone();
1279  finalpolys->push_back(newGeom);
1280  }
1281 
1282  p_polygons = globalFactory.createMultiPolygon(*finalpolys);
1283 
1284  delete finalpolys;
1285  delete newGeom;
1286  delete newLonLatPts;
1287  delete pts;
1288  delete pts2;
1289  }
1290  catch(geos::util::IllegalArgumentException *geosIll) {
1291  std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1292  msg += "geos illegal argument [" + IString(geosIll->what()) + "]";
1293  delete geosIll;
1295  }
1296  catch(geos::util::GEOSException *geosExc) {
1297  std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1298  msg += "geos exception [" + IString(geosExc->what()) + "]";
1299  delete geosExc;
1301  }
1302  catch(IException &e) {
1303  std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1304  msg += "isis operation exception [" + IString(e.what()) + "]";
1305  throw IException(e, IException::Unknown, msg, _FILEINFO_);
1306  }
1307  catch(...) {
1308  std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1309  msg += "unknown exception";
1311  }
1312  }
1313 
1314 
1323  void ImagePolygon::ReadData(std::istream &is) {
1324 
1325  streampos sbyte = p_startByte - 1;
1326  is.seekg(sbyte, std::ios::beg);
1327  if (!is.good()) {
1328  QString msg = "Error preparing to read data from " + p_type +
1329  " [" + p_blobName + "]";
1330  throw IException(IException::Io, msg, _FILEINFO_);
1331  }
1332 
1333  char *buf = new char[p_nbytes+1];
1334  memset(buf, 0, p_nbytes + 1);
1335 
1336  is.read(buf, p_nbytes);
1337 
1338  p_polyStr = buf;
1339 
1340  delete [] buf;
1341 
1342  if (!is.good()) {
1343  QString msg = "Error reading data from " + p_type + " [" +
1344  p_blobName + "]";
1345  throw IException(IException::Io, msg, _FILEINFO_);
1346  }
1347 
1348  geos::io::WKTReader *wkt = new geos::io::WKTReader(&(globalFactory));
1350  delete wkt;
1351  }
1352 
1353 
1356  geos::io::WKTWriter *wkt = new geos::io::WKTWriter();
1357 
1358  // Check to see p_polygons is valid data
1359  if (!p_polygons) {
1360  string msg = "Cannot write a NULL polygon!";
1362  }
1363  p_polyStr = wkt->write(p_polygons);
1364  p_nbytes = p_polyStr.size();
1365 
1366  delete wkt;
1367  }
1368 
1369 
1377  void ImagePolygon::WriteData(std::fstream &os) {
1378  os.write(p_polyStr.c_str(), p_nbytes);
1379  }
1380 
1381 
1390  double ImagePolygon::DistanceSquared(const geos::geom::Coordinate *p1, const geos::geom::Coordinate *p2) {
1391  return ((p2->x - p1->x) * (p2->x - p1->x)) + ((p2->y - p1->y) * (p2->y - p1->y));
1392  }
1393 
1394 
1401  bool hasFourCorners = true;
1402  hasFourCorners &= SetImage(1, 1);
1403  hasFourCorners &= SetImage(p_cubeSamps, 1);
1404  hasFourCorners &= SetImage(p_cubeSamps, p_cubeLines);
1405  hasFourCorners &= SetImage(1, p_cubeLines);
1406  return !hasFourCorners;
1407  }
1408 
1409 
1410 
1424  geos::geom::Coordinate ImagePolygon::FindBestPoint(geos::geom::Coordinate *currentPoint,
1425  geos::geom::Coordinate newPoint,
1426  geos::geom::Coordinate lastPoint) {
1427  geos::geom::Coordinate result = newPoint;
1428  // Use a binary search to snap to the limb when needed
1429  if (p_sampinc > 1 || p_lineinc > 1) {
1430 
1431  // Pull the invalid point back inside the image
1432  double x = lastPoint.x;
1433  double y = lastPoint.y;
1434  if (x < p_cubeStartSamp) {
1435  x = p_cubeStartSamp;
1436  }
1437  else if (x > p_cubeSamps) {
1438  x = p_cubeSamps;
1439  }
1440  if (y < p_cubeStartLine) {
1441  y = p_cubeStartLine;
1442  }
1443  else if (y > p_cubeLines) {
1444  y = p_cubeLines;
1445  }
1446  geos::geom::Coordinate invalid(x, y);
1447  geos::geom::Coordinate valid(result.x, result.y);
1448 
1449  // Find the best valid Coordinate
1450  while (!SetImage(invalid.x, invalid.y)) {
1451  int x, y;
1452  if (invalid.x > valid.x) {
1453  x = (int)invalid.x - 1;
1454  }
1455  else if (invalid.x < valid.x) {
1456  x = (int)invalid.x + 1;
1457  }
1458  else {
1459  x = (int)invalid.x;
1460  }
1461  if (invalid.y > valid.y) {
1462  y = (int)invalid.y - 1;
1463  }
1464  else if (invalid.y < valid.y) {
1465  y = (int)invalid.y + 1;
1466  }
1467  else {
1468  y = (int)invalid.y;
1469  }
1470  invalid = geos::geom::Coordinate(x, y);
1471  }
1472 
1473  result = FixCornerSkip(currentPoint, invalid);
1474  }
1475 
1476  return result;
1477  }
1478 
1479 
1489  geos::geom::Coordinate ImagePolygon::FixCornerSkip(geos::geom::Coordinate *currentPoint,
1490  geos::geom::Coordinate newPoint) {
1491  geos::geom::Coordinate originalPoint = newPoint;
1492  geos::geom::Coordinate modPoint = newPoint;
1493 
1494  // Step Too big
1496  return newPoint;
1497  }
1498 
1499  // An upper left corner
1500  else if (currentPoint->x < newPoint.x && currentPoint->y > newPoint.y) {
1501  while (newPoint.x >= currentPoint->x && SetImage(newPoint.x, newPoint.y)) {
1502  modPoint = newPoint;
1503  newPoint.x -= 1;
1504  }
1505  }
1506 
1507  // An upper right corner
1508  else if (currentPoint->y < newPoint.y && currentPoint->x < newPoint.x) {
1509  while (newPoint.y >= currentPoint->y && SetImage(newPoint.x, newPoint.y)) {
1510  modPoint = newPoint;
1511  newPoint.y -= 1;
1512  }
1513  }
1514 
1515  // An lower right corner
1516  else if (currentPoint->x > newPoint.x && currentPoint->y < newPoint.y) {
1517  while (newPoint.x <= currentPoint->x && SetImage(newPoint.x, newPoint.y)) {
1518  modPoint = newPoint;
1519  newPoint.x += 1;
1520  }
1521  }
1522 
1523  // An lower left corner
1524  else if (currentPoint->y > newPoint.y && currentPoint->x > newPoint.x) {
1525  while (newPoint.y <= currentPoint->y && SetImage(newPoint.x, newPoint.y)) {
1526  modPoint = newPoint;
1527  newPoint.y += 1;
1528  }
1529  }
1530 
1531  if (currentPoint->x == modPoint.x && currentPoint->y == modPoint.y) {
1532  return originalPoint;
1533  }
1534  return modPoint;
1535  }
1536 
1537 
1545  void ImagePolygon::FindSubpixel(std::vector<geos::geom::Coordinate> & points) {
1546  if (p_subpixelAccuracy > 0) {
1547 
1548  // Fix the polygon with subpixel accuracy
1549  geos::geom::Coordinate old = points.at(0);
1550  bool didStartingPoint = false;
1551  for (unsigned int pt = 1; !didStartingPoint; pt ++) {
1552  if (pt >= points.size() - 1) {
1553  pt = 0;
1554  didStartingPoint = true;
1555  }
1556 
1557  // Binary Coordinate Search
1558  double maxStep = std::max(p_sampinc, p_lineinc);
1559  double stepY = (old.x - points.at(pt + 1).x) / maxStep;
1560  double stepX = (points.at(pt + 1).y - old.y) / maxStep;
1561 
1562  geos::geom::Coordinate valid = points.at(pt);
1563  geos::geom::Coordinate invalid(valid.x + stepX, valid.y + stepY);
1564 
1565  for (int itt = 0; itt < p_subpixelAccuracy; itt ++) {
1566  geos::geom::Coordinate half((valid.x + invalid.x) / 2.0, (valid.y + invalid.y) / 2.0);
1567  if (SetImage(half.x, half.y) && InsideImage(half.x, half.y)) {
1568  valid = half;
1569  }
1570  else {
1571  invalid = half;
1572  }
1573  }
1574 
1575  old = points.at(pt);
1576 
1577  // Set new coordinate
1578  points[pt] = valid;
1579 
1580  }
1581 
1582  // Fix starting point
1583  points[points.size()-1] = geos::geom::Coordinate(points[0].x, points[0].y);
1584 
1585  }
1586  }
1587 
1588 
1589 } // end namespace isis
1590 
double validSampleDim()
Retuns the maximum valid sample width of the cube set with either initCube() or Create() ...
const char * what() const
Returns a string representation of this exception in its current state.
Definition: IException.cpp:391
QString p_type
Type of data stored in the buffer.
Definition: Blob.h:100
geos::geom::Coordinate FindFirstPoint()
Finds the first point that projects in an image.
void IgnoreElevationModel(bool ignore)
This allows you to ignore the cube elevation model and use the ellipse.
Definition: Sensor.cpp:75
const double Null
Value for an Isis Null pixel.
Definition: SpecialPixel.h:110
double validLineDim()
Retuns the maximum valid line width of the cube set with either initCube() or Create() ...
geos::geom::Coordinate FixCornerSkip(geos::geom::Coordinate *currentPoint, geos::geom::Coordinate newPoint)
Looks at the next possible point relative to the lasts and attempts to adjust the point outward to gr...
geos::geom::Coordinate * m_rightCoord
The cube&#39;s right-most valid coord.
Definition: ImagePolygon.h:286
geos::geom::Coordinate * m_topCoord
The cube&#39;s top-most valid coord.
Definition: ImagePolygon.h:287
bool InsideImage(double sample, double line)
This returns true if sample/line are inside the cube.
int p_subpixelAccuracy
The subpixel accuracy to use.
Definition: ImagePolygon.h:302
Universal Ground Map.
QString p_blobName
Name of the Blob object.
Definition: Blob.h:95
Camera * camera()
Return a camera associated with the cube.
Definition: Cube.cpp:1166
int p_cubeLines
The number of lines in the cube.
Definition: ImagePolygon.h:293
int sampleCount() const
Definition: Cube.cpp:1452
bool HasProjection()
Checks to see if the camera object has a projection.
Definition: Camera.cpp:2650
Buffer for containing a three dimensional section of an image.
Definition: Brick.h:61
bool p_ellipsoid
Uses an ellipsoid if a limb is detected.
Definition: ImagePolygon.h:300
void append(const IException &exceptionSource)
Appends the given exception (and its list of previous exceptions) to this exception&#39;s causational exc...
Definition: IException.cpp:425
Namespace for the standard library.
void SetBasePosition(const int start_sample, const int start_line, const int start_band)
This method is used to set the base position of the shape buffer.
Definition: Brick.h:136
int p_nbytes
Size of blob data (in bytes)
Definition: Blob.h:99
Camera * initCube(Cube &cube, int ss=1, int sl=1, int ns=0, int nl=0, int band=1)
Create a Polygon from given cube.
bool p_isProjected
True when the provided cube is projected.
Definition: ImagePolygon.h:273
Projection * projection()
Definition: Cube.cpp:1439
geos::geom::Coordinate * m_leftCoord
The cube&#39;s left-most valid coord.
Definition: ImagePolygon.h:285
This error is for when a programmer made an API call that was illegal.
Definition: IException.h:162
A type of error that occurred when performing an actual I/O operation.
Definition: IException.h:171
void Create(Cube &cube, int sinc=1, int linc=1, int ss=1, int sl=1, int ns=0, int nl=0, int band=1, bool increasePrecision=false)
Create a Polygon from given cube.
double p_emission
The maximum emission angle to consider valid.
Definition: ImagePolygon.h:298
Cube * p_cube
The cube provided.
Definition: ImagePolygon.h:272
void WalkPoly()
Walks the image finding its lon lat polygon and stores it to p_pts.
static geos::geom::Geometry * Intersect(const geos::geom::Geometry *geom1, const geos::geom::Geometry *geom2)
This applies the geos Intersect operator.
UniversalGroundMap * p_gMap
The cube&#39;s ground map.
Definition: ImagePolygon.h:283
double EmissionAngle() const
Returns the emission angle in degrees.
Definition: Sensor.cpp:339
double WorldX() const
This returns the world X coordinate provided SetGround, SetCoordinate, SetUniversalGround, or SetWorld returned with success.
Definition: Projection.cpp:539
Brick * p_brick
Used to check for valid DNs.
Definition: ImagePolygon.h:275
geos::geom::Coordinate * m_botCoord
The cube&#39;s bot-most valid coord.
Definition: ImagePolygon.h:288
geos::geom::Coordinate FindNextPoint(geos::geom::Coordinate *currentPoint, geos::geom::Coordinate lastPoint, int recursionDepth=0)
Finds the next point on the image using a left hand rule walking algorithm.
double DistanceSquared(const geos::geom::Coordinate *p1, const geos::geom::Coordinate *p2)
Calculates the distance squared between two coordinates.
int p_lineinc
The increment for walking along the polygon in the line direcction.
Definition: ImagePolygon.h:296
std::string p_polyStr
The string representation of the polygon.
Definition: ImagePolygon.h:281
void WriteInit()
Initializes for writing polygon to cube blob.
static geos::geom::MultiPolygon * MakeMultiPolygon(const geos::geom::Geometry *geom)
Make a geos::geom::MultiPolygon out of the components of the argument.
#define _FILEINFO_
Macro for the filename and line number.
Definition: IException.h:40
bool SetImage(const double sample, const double line)
Sets the sample/line values of the cube to get lat/lon values.
double Sample()
Returns the current sample number.
Definition: Camera.cpp:2702
A type of error that could only have occurred due to a mistake on the user&#39;s part (e...
Definition: IException.h:142
void ReadData(std::istream &is)
Reads Multipolygon from cube blob.
A type of error that cannot be classified as any of the other error types.
Definition: IException.h:134
Isis::Projection * Projection() const
Return the projection associated with the ground map (NULL implies none)
double UniversalLatitude() const
Returns the universal latitude of the camera model or projection.
void read(Blob &blob) const
This method will read data from the specified Blob object.
Definition: Cube.cpp:724
int p_cubeStartLine
The the line of the first valid point in the cube.
Definition: ImagePolygon.h:291
geos::geom::MultiPolygon * p_polygons
The multipolygon of the image.
Definition: ImagePolygon.h:279
double p_incidence
The maximum incidence angle to consider valid.
Definition: ImagePolygon.h:299
int p_cubeStartSamp
The the sample of the first valid point in the cube.
Definition: ImagePolygon.h:290
void BasicMapping(Pvl &map)
Writes the basic mapping group to the specified Pvl.
Definition: Camera.cpp:1375
void FixPolePoly(std::vector< geos::geom::Coordinate > *crossingPoints)
If the cube crosses the 0/360 boundary and contains a pole, Some points are added to allow the polygo...
Container for cube-like labels.
Definition: Pvl.h:135
PixelType pixelType() const
Definition: Cube.cpp:1403
bool IsLimb()
Returns True when the input image is a limb image.
int p_cubeSamps
The number of samples in the cube.
Definition: ImagePolygon.h:292
bool SetImage(double sample, double line)
Returns whether the sample/line postion was set successfully in the camera model or projection...
void print() const
Prints a string representation of this exception to stderr.
Definition: IException.cpp:461
void SetBand(const int band)
Set the image band number.
void Fix360Poly()
If the cube crosses the 0/360 boundary and does not include a pole, the polygon is separated into mul...
virtual QString fileName() const
Returns the opened cube&#39;s filename.
Definition: Cube.cpp:1208
double Line()
Returns the current line number.
Definition: Camera.cpp:2722
int lineCount() const
Definition: Cube.cpp:1379
void WriteData(std::fstream &os)
Writes polygon to cube blob.
void calcImageBorderCoordinates()
Calculates the four border points, particularly for validSampleDim() and validLineDim() ...
Isis exception class.
Definition: IException.h:107
Adds specific functionality to C++ strings.
Definition: IString.h:181
void MoveBackInsideImage(double &sample, double &line, double sinc, double linc)
This method ensures sample/line after sinc/linc have been applied is inside the image.
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31
Isis::Camera * Camera() const
Return the camera associated with the ground map (NULL implies none)
bool IsNullPixel(const double d)
Returns if the input pixel is null.
Definition: SpecialPixel.h:250
geos::geom::Coordinate FindBestPoint(geos::geom::Coordinate *currentPoint, geos::geom::Coordinate newPoint, geos::geom::Coordinate lastPoint)
While walking the image in sample/line space, this function finds the best valid point between the fi...
bool SetUniversalGround(double lat, double lon)
Returns whether the lat/lon position was set successfully in the camera model or projection.
double IncidenceAngle() const
Returns the incidence angle in degrees.
Definition: Sensor.cpp:350
geos::geom::CoordinateSequence * p_pts
The sequence of coordinates that compose the boundry of the image.
Definition: ImagePolygon.h:277
void FindSubpixel(std::vector< geos::geom::Coordinate > &points)
Takes p_polygons in sample/line space and finds its subpixel accuracy.
double UniversalLongitude() const
Returns the universal longitude of the camera model or projection.
BigInt p_startByte
Byte blob data starts at in buffer.
Definition: Blob.h:98
~ImagePolygon()
Destroys the Polygon object.
double WorldY() const
This returns the world Y coordinate provided SetGround, SetCoordinate, SetUniversalGround, or SetWorld returned with success.
Definition: Projection.cpp:559
int p_sampinc
The increment for walking along the polygon in the sample direction.
Definition: ImagePolygon.h:295
IO Handler for Isis Cubes.
Definition: Cube.h:170