Isis 3 Programmer Reference
ImagePolygon.cpp
1
6/* SPDX-License-Identifier: CC0-1.0 */
7
8#include <string>
9#include <iostream>
10#include <vector>
11
12#include <QDebug>
13
14#include <geos/geom/Geometry.h>
15#include <geos/geom/Polygon.h>
16#include <geos/geom/CoordinateSequence.h>
17#include <geos/algorithm/LineIntersector.h>
18#include <geos/util/IllegalArgumentException.h>
19#include <geos/util/TopologyException.h>
20#include <geos/util/GEOSException.h>
21#include <geos/io/WKTReader.h>
22#include <geos/io/WKTWriter.h>
23#include <geos/operation/distance/DistanceOp.h>
24
25#include "ImagePolygon.h"
26#include "IString.h"
27#include "SpecialPixel.h"
28#include "PolygonTools.h"
29
30using namespace std;
31
32namespace Isis {
33
39 p_polygons = NULL;
40
41 p_cube = NULL;
42
43 m_leftCoord = NULL;
44 m_rightCoord = NULL;
45 m_topCoord = NULL;
46 m_botCoord = NULL;
47
50
51 p_emission = 180.0;
52 p_incidence = 180.0;
53
54 p_subpixelAccuracy = 50; //An accuracte and quick number
55
56 p_ellipsoid = false;
57 }
58
59
65 p_polyStr = string(blob.getBuffer(), blob.Size());
66
67 geos::io::WKTReader *wkt = new geos::io::WKTReader(&(*globalFactory));
69
70 p_pts = new geos::geom::CoordinateSequence;
71
72// for (auto poly : *p_polygons) {
73 for(unsigned int g = 0; g < p_polygons->getNumGeometries(); ++g) {
74 const geos::geom::Polygon *poly =
75 dynamic_cast<const geos::geom::Polygon *>(p_polygons->getGeometryN(g));
76 geos::geom::CoordinateSequence coordArray = geos::geom::CoordinateSequence(*(poly->getCoordinates()));
77 for (size_t i = 0; i < coordArray.getSize(); i++) {
78 p_pts->add(geos::geom::Coordinate(coordArray.getAt(i)));
79 }
80 }
81 }
82
83
86 delete p_polygons;
87 p_polygons = NULL;
88
89 p_cube = NULL;
90
91 delete m_leftCoord;
92 m_leftCoord = NULL;
93
94 delete m_rightCoord;
95 m_rightCoord = NULL;
96
97 delete m_topCoord;
98 m_topCoord = NULL;
99
100 delete m_botCoord;
101 m_botCoord = NULL;
102 }
103
104
120 Camera * ImagePolygon::initCube(Cube &cube, int ss, int sl,
121 int ns, int nl, int band) {
122 p_gMap = new UniversalGroundMap(cube);
123 p_gMap->SetBand(band);
124
125 p_cube = &cube;
126
127 Camera *cam = NULL;
128 p_isProjected = false;
129
130 try {
131 cam = cube.camera();
132 }
133 catch(IException &camError) {
134 try {
135 cube.projection();
136 p_isProjected = true;
137 }
138 catch(IException &projError) {
139 QString msg = "Can not create polygon, ";
140 msg += "cube [" + cube.fileName();
141 msg += "] is not a camera or map projection";
142
143 IException polyError(IException::User, msg, _FILEINFO_);
144
145 polyError.append(camError);
146 polyError.append(projError);
147
148 throw polyError;
149 }
150 }
151 if (cam != NULL) p_isProjected = cam->HasProjection();
152
153 // Create brick for use in SetImage
154 p_brick = new Brick(1, 1, 1, cube.pixelType());
155
156 //------------------------------------------------------------------------
157 // Save cube number of samples and lines for later use.
158 //------------------------------------------------------------------------
159 p_cubeSamps = cube.sampleCount();
160 p_cubeLines = cube.lineCount();
161
162 if (ns != 0) {
163 p_cubeSamps = std::min(p_cubeSamps, ss + ns);
164 }
165
166 if (nl != 0) {
167 p_cubeLines = std::min(p_cubeLines, sl + nl);
168 }
169
170 p_cubeStartSamp = ss;
171 p_cubeStartLine = sl;
172
173 if (p_ellipsoid && IsLimb() && p_gMap->Camera()) {
174 try {
176 }
177 catch(IException &) {
178 std::string msg = "Cannot use an ellipsoid shape model";
179 msg += " on a limb image without a camera.";
180 throw IException(IException::User, msg, _FILEINFO_);
181 }
182 }
183
184 return cam;
185 }
186
187
216 void ImagePolygon::Create(Cube &cube, int sinc, int linc,
217 int ss, int sl, int ns, int nl, int band,
218 bool increasePrecision) {
219
220 Camera *cam = NULL;
221
222 if (sinc < 1 || linc < 1) {
223 std::string msg = "Sample and line increments must be 1 or greater";
224 throw IException(IException::User, msg, _FILEINFO_);
225 }
226
227 cam = initCube(cube, ss, sl, ns, nl, band);
228
229 // Reduce the increment size to find a valid polygon
230 bool polygonGenerated = false;
231 while (!polygonGenerated) {
232 try {
233 p_sampinc = sinc;
234 p_lineinc = linc;
235
236 p_pts = NULL;
237 p_pts = new geos::geom::CoordinateSequence();
238
239 WalkPoly();
240
241 polygonGenerated = true;
242 }
243 catch (IException &e) {
244 delete p_pts;
245
246 sinc = sinc * 2 / 3;
247 linc = linc * 2 / 3;
248
249 if (increasePrecision && (sinc > 1 || linc > 1)) {
250 }
251 else {
252 e.print(); // This should be a NAIF error
253 QString msg = "Cannot find polygon for image "
254 "[" + cube.fileName() + "]: ";
255 msg += increasePrecision ? "Cannot increase precision any further" :
256 "The increment/step size might be too large";
257 throw IException(IException::User, msg, _FILEINFO_);
258 }
259 }
260 }
261
262 /*------------------------------------------------------------------------
263 / If image contains 0/360 boundary, the polygon needs to be split up
264 / into multi polygons.
265 /-----------------------------------------------------------------------*/
266 if (cam) {
267 Pvl defaultMap;
268 cam->BasicMapping(defaultMap);
269 }
270
271 // Create the polygon, fixing if needed
272 Fix360Poly();
273
274 if (p_brick != 0) delete p_brick;
275
276 if (p_gMap->Camera())
278 }
279
280
281 void ImagePolygon::Create(std::vector<std::vector<double>> polyCoordinates) {
282 p_pts = new geos::geom::CoordinateSequence();
283
284 for (std::vector<double> coord : polyCoordinates) {
285 p_pts->add(geos::geom::Coordinate(coord[0], coord[1]));
286 }
287
288 std::vector<const geos::geom::Geometry *> *polys = new std::vector<const geos::geom::Geometry *>;
289 const geos::geom::Polygon *poly = globalFactory->createPolygon(globalFactory->createLinearRing(*p_pts)).release();
290 polys->push_back(poly->clone().release());
291 // createMultiPolygon takes ownership of polys argument
292 p_polygons = globalFactory->createMultiPolygon(*polys).release();
293
294 Fix360Poly();
295 }
296
297
314 geos::geom::Coordinate ImagePolygon::FindNextPoint(geos::geom::Coordinate *currentPoint,
315 geos::geom::Coordinate lastPoint,
316 int recursionDepth) {
317 double x = lastPoint.x - currentPoint->x;
318 double y = lastPoint.y - currentPoint->y;
319 geos::geom::Coordinate result;
320
321 // Check to see if depth is too deep (walked all the way around and found nothing)
322 if (recursionDepth > 6) {
323 return *currentPoint;
324 }
325
326 // Check and walk in appropriate direction
327 if (x == 0.0 && y == 0.0) { // Find the starting point on an image
328 for (int line = -1 * p_lineinc; line <= 1 * p_lineinc; line += p_lineinc) {
329 for (int samp = -1 * p_sampinc; samp <= 1 * p_sampinc; samp += p_sampinc) {
330 double s = currentPoint->x + samp;
331 double l = currentPoint->y + line;
332 // Try the next left hand rule point if (s,l) does not produce a
333 // lat/long or it is not on the image.
334 if (!InsideImage(s, l) || !SetImage(s, l)) {
335 geos::geom::Coordinate next(s, l);
336 return FindNextPoint(currentPoint, next);
337 }
338 }
339 }
340
341 std::string msg = "Unable to create image footprint. Starting point is not on the edge of the image.";
342 throw IException(IException::Programmer, msg, _FILEINFO_);
343 }
344 else if (x < 0 && y < 0) { // current is top left
345 geos::geom::Coordinate next(currentPoint->x, currentPoint->y - 1 * p_lineinc);
346 MoveBackInsideImage(next.x, next.y, 0, -p_lineinc);
347 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
348 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
349 }
350 else {
351 result = FindBestPoint(currentPoint, next, lastPoint);
352 }
353 }
354 else if (x == 0.0 && y < 0) { // current is top
355 geos::geom::Coordinate next(currentPoint->x + 1 * p_sampinc, currentPoint->y - 1 * p_lineinc);
356 MoveBackInsideImage(next.x, next.y, p_sampinc, -p_lineinc);
357 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
358 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
359 }
360 else {
361 result = FindBestPoint(currentPoint, next, lastPoint);
362 }
363 }
364 else if (x > 0 && y < 0) { // current is top right
365 geos::geom::Coordinate next(currentPoint->x + 1 * p_sampinc, currentPoint->y);
366 MoveBackInsideImage(next.x, next.y, p_sampinc, 0);
367 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
368 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
369 }
370 else {
371 result = FindBestPoint(currentPoint, next, lastPoint);
372 }
373 }
374 else if (x > 0 && y == 0.0) { // current is right
375 geos::geom::Coordinate next(currentPoint->x + 1 * p_sampinc, currentPoint->y + 1 * p_lineinc);
376 MoveBackInsideImage(next.x, next.y, p_sampinc, p_lineinc);
377 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
378 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
379 }
380 else {
381 result = FindBestPoint(currentPoint, next, lastPoint);
382 }
383 }
384 else if (x > 0 && y > 0) { // current is bottom right
385 geos::geom::Coordinate next(currentPoint->x, currentPoint->y + 1 * p_lineinc);
386 MoveBackInsideImage(next.x, next.y, 0, p_lineinc);
387 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
388 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
389 }
390 else {
391 result = FindBestPoint(currentPoint, next, lastPoint);
392 }
393 }
394 else if (x == 0.0 && y > 0) { // current is bottom
395 geos::geom::Coordinate next(currentPoint->x - 1 * p_sampinc, currentPoint->y + 1 * p_lineinc);
396 MoveBackInsideImage(next.x, next.y, -p_sampinc, p_lineinc);
397 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
398 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
399 }
400 else {
401 result = FindBestPoint(currentPoint, next, lastPoint);
402 }
403 }
404 else if (x < 0 && y > 0) { // current is bottom left
405 geos::geom::Coordinate next(currentPoint->x - 1 * p_sampinc, currentPoint->y);
406 MoveBackInsideImage(next.x, next.y, -p_sampinc, 0);
407 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
408 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
409 }
410 else {
411 result = FindBestPoint(currentPoint, next, lastPoint);
412 }
413 }
414 else if (x < 0 && y == 0.0) { // current is left
415 geos::geom::Coordinate next(currentPoint->x - 1 * p_sampinc, currentPoint->y - 1 * p_lineinc);
416 MoveBackInsideImage(next.x, next.y, -p_sampinc, -p_lineinc);
417 if (!recursionDepth || !InsideImage(next.x, next.y) || !SetImage(next.x , next.y)) {
418 result = FindNextPoint(currentPoint, next, recursionDepth + 1);
419 }
420 else {
421 result = FindBestPoint(currentPoint, next, lastPoint);
422 }
423 }
424 else {
425 std::string msg = "Unable to create image footprint. Error walking image.";
426 throw IException(IException::Unknown, msg, _FILEINFO_);
427 }
428
429
430 return result;
431 }
432
433
444 void ImagePolygon::MoveBackInsideImage(double &sample, double &line, double sinc, double linc) {
445 // snap to centers of pixels!
446
447 // Starting sample to snap to
448 const double startSample = p_cubeStartSamp;
449 // Ending sample to snap to
450 const double endSample = p_cubeSamps;
451 // Starting line to snap to
452 const double startLine = p_cubeStartLine;
453 // Ending line to snap to
454 const double endLine = p_cubeLines;
455 // Original sample for this point (before sinc)
456 const double origSample = sample - sinc;
457 // Original line for this point (before linc)
458 const double origLine = line - linc;
459
460 // We moved left off the image - snap to the edge
461 if (sample < startSample && sinc < 0) {
462 // don't snap if we started at the edge
463 if (origSample == startSample) {
464 return;
465 }
466
467 sample = startSample;
468 }
469
470 // We moved right off the image - snap to the edge
471 if (sample > endSample && sinc > 0) {
472 // don't snap if we started at the edge
473 if (origSample == endSample) {
474 return;
475 }
476
477 sample = endSample;
478 }
479
480 // We moved up off the image - snap to the edge
481 if (line < startLine && linc < 0) {
482 // don't snap if we started at the edge
483 if (origLine == startLine) {
484 return;
485 }
486
487 line = startLine;
488 }
489
490 // We moved down off the image - snap to the edge
491 if (line > endLine && linc > 0) {
492 // don't snap if we started at the edge
493 if (fabs(origLine - endLine) < 0.5) {
494 return;
495 }
496
497 line = endLine;
498 }
499
500 return;
501 }
502
503
512 bool ImagePolygon::InsideImage(double sample, double line) {
513 return (sample >= p_cubeStartSamp - 0.5 &&
514 line > p_cubeStartLine - 0.5 &&
515 sample <= p_cubeSamps + 0.5 &&
516 line <= p_cubeLines + 0.5);
517 }
518
519
526 geos::geom::Coordinate ImagePolygon::FindFirstPoint() {
527 // @todo: Brute force method, should be improved
528 for (int sample = p_cubeStartSamp; sample <= p_cubeSamps; sample++) {
529 for (int line = p_cubeStartLine; line <= p_cubeLines; line++) {
530 if (SetImage(sample, line)) {
531 // An outlier check. Make sure that the pixel we use to start
532 // constructing a polygon is not surrounded by a bunch of invalid
533 // positions.
534 geos::geom::Coordinate firstPoint(sample, line);
535 geos::geom::Coordinate lastPoint = firstPoint;
536 if (!firstPoint.equals(FindNextPoint(&firstPoint, lastPoint))) {
537 return firstPoint;
538 }
539 }
540 }
541 }
542
543 // Check to make sure a point was found
544 std::string msg = "No lat/lon data found for image";
545 throw IException(IException::User, msg, _FILEINFO_);
546 }
547
548
554 double result = 0.0;
555
558 result = m_rightCoord->x - m_leftCoord->x + 1;
559
560 return result;
561 }
562
563
569 double result = 0.0;
570
572 if (m_topCoord && m_botCoord)
573 result = m_botCoord->y - m_topCoord->y + 1;
574
575 return result;
576 }
577
578
588
589 for (int line = p_cubeStartLine; !m_leftCoord && line <= p_cubeLines; line++) {
590 for (int sample = p_cubeStartSamp; !m_leftCoord && sample <= p_cubeSamps; sample++) {
591 if (SetImage(sample, line)) {
592 m_leftCoord = new geos::geom::Coordinate(sample, line);
593 }
594 }
595 }
596
597 if (m_leftCoord) {
598 for (int line = p_cubeStartLine; !m_rightCoord && line <= p_cubeLines; line++) {
599 for (int sample = p_cubeSamps; !m_rightCoord && sample >= m_leftCoord->x; sample--) {
600 if (SetImage(sample, line)) {
601 m_rightCoord = new geos::geom::Coordinate(sample, line);
602 }
603 }
604 }
605 }
606
607 if (m_leftCoord && m_rightCoord) {
608 for (int sample = (int)m_leftCoord->x; !m_topCoord && sample <= m_rightCoord->x; sample++) {
609 for (int line = 1; !m_topCoord && line <= p_cubeLines; line++) {
610 if (SetImage(sample, line)) {
611 m_topCoord = new geos::geom::Coordinate(sample, line);
612 }
613 }
614 }
615 }
616
618 for (int sample = (int)m_leftCoord->x; !m_botCoord && sample <= m_rightCoord->x; sample++) {
619 for (int line = p_cube->lineCount(); !m_botCoord && line >= m_topCoord->y; line--) {
620 if (SetImage(sample, line)) {
621 m_botCoord = new geos::geom::Coordinate(sample, line);
622 }
623 }
624 }
625 }
626 }
627
628
636 vector<geos::geom::Coordinate> points;
637 double lat, lon, prevLat, prevLon;
638
639 // Find the edge of the polygon
640 geos::geom::Coordinate firstPoint = FindFirstPoint();
641
642 points.push_back(firstPoint);
643 //************************
644 // Start walking the edge
645 //************************
646
647 // set up for intialization
648 geos::geom::Coordinate currentPoint = firstPoint;
649 geos::geom::Coordinate lastPoint = firstPoint;
650 geos::geom::Coordinate tempPoint;
651
652 do {
653 tempPoint = FindNextPoint(&currentPoint, lastPoint);
654
655 // First check if the distance is within range of skipping
656 bool snapToFirstPoint = true;
657
658 // Never needs to find firstPoint on incements of 1
659 snapToFirstPoint &= (p_sampinc != 1 && p_lineinc != 1);
660
661 // Prevents catching the first point as the last
662 snapToFirstPoint &= (points.size() > 2);
663
664 // This method fails for steps larger than line/sample length
665 snapToFirstPoint &= (p_sampinc < p_cubeSamps && p_lineinc < p_cubeLines);
666
667 // Checks for appropriate distance without sqrt() call
668 int minStepSize = std::min(p_sampinc, p_lineinc);
669 snapToFirstPoint &= (DistanceSquared(&currentPoint, &firstPoint) < (minStepSize * minStepSize));
670
671 // Searches for skipped firstPoint due to a large pixinc, by using a pixinc of 1
672 if (snapToFirstPoint) {
673 tempPoint = firstPoint;
674 }
675 // If pixinc is greater than line or sample dimention, then we could
676 // skip firstPoint. This checks for that case.
677 else if (p_sampinc > p_cubeSamps || p_lineinc > p_cubeLines) {
678 // This is not expensive because incement must be large
679 for (int pt = 0; pt < (int)points.size(); pt ++) {
680 if (points[pt].equals(tempPoint)) {
681 tempPoint = firstPoint;
682 break;
683 }
684 }
685 }
686 // Failed to find the next point
687 if (tempPoint.equals(currentPoint)) {
688 geos::geom::Coordinate oldDuplicatePoint = tempPoint;
689
690 // Init vars for first run through the loop
691 tempPoint = lastPoint;
692 lastPoint = currentPoint;
693 currentPoint = tempPoint;
694
695 // Must be 3 (not 2) to prevent the revisit of the starting point,
696 // resulting in an infinite loop
697 if (points.size() < 3) {
698 std::string msg = "Failed to find next point in the image.";
699 throw IException(IException::Programmer, msg, _FILEINFO_);
700 }
701
702 // remove last point from the list
703 points.pop_back();
704
705 tempPoint = FindNextPoint(&currentPoint, lastPoint, 1);
706
707 if (tempPoint.equals(currentPoint) || tempPoint.equals(oldDuplicatePoint)) {
708 std::string msg = "Failed to find next valid point in the image.";
709 throw IException(IException::Programmer, msg, _FILEINFO_);
710 }
711 }
712
713
714 // Check for triangle cycles and try to fix
715 if ((p_sampinc > 1 || p_lineinc > 1) && points.size() >= 3) {
716 if (points[points.size()-3].x == tempPoint.x &&
717 points[points.size()-3].y == tempPoint.y) {
718 // Remove the triangle from the list
719 points.pop_back();
720 points.pop_back();
721 points.pop_back();
722 // Reset the current (to be last) point
723 currentPoint = points[points.size()-1];
724 // change increment to prevent randomly bad pixels in the image
725 if (p_sampinc > 1) p_sampinc --;
726 if (p_lineinc > 1) p_lineinc --;
727 }
728
735 if (points.size() > 250) {
736 int cycleStart = 0;
737 int cycleEnd = 0;
738
739 for (unsigned int pt = 1; pt < points.size() && cycleStart == 0; pt ++) {
740 for (unsigned int check = pt + 1; check < points.size() && cycleStart == 0; check ++) {
741 if (points[pt] == points[check]) {
742 cycleStart = pt;
743 cycleEnd = check;
744 }
745 }
746 }
747
748 // If a cycle was found, make it the polygon
749 if (cycleStart != 0) {
750 vector<geos::geom::Coordinate> cyclePoints;
751 for (int pt = cycleStart; pt <= cycleEnd; pt ++) {
752 cyclePoints.push_back(points[pt]);
753 }
754
755 points = cyclePoints;
756 break;
757 }
758 }
759
760 }
761
762 lastPoint = currentPoint;
763 currentPoint = tempPoint;
764 points.push_back(currentPoint);
765
766 }
767 while (!currentPoint.equals(firstPoint));
768
769 if (points.size() <= 3) {
770 std::string msg = "Failed to find enough points on the image.";
771 throw IException(IException::Programmer, msg, _FILEINFO_);
772 }
773
774 FindSubpixel(points);
775
776 prevLat = 0;
777 prevLon = 0;
778 // this vector stores crossing points, where the image crosses the
779 // meridian. It stores the first coordinate of the pair in its vector
780 vector<geos::geom::Coordinate> *crossingPoints = new vector<geos::geom::Coordinate>;
781 for (unsigned int i = 0; i < points.size(); i++) {
782 geos::geom::Coordinate *temp = &(points.at(i));
783 SetImage(temp->x, temp->y);
784 lon = p_gMap->UniversalLongitude();
785 lat = p_gMap->UniversalLatitude();
786 if (abs(lon - prevLon) >= 180 && i != 0) {
787 crossingPoints->push_back(geos::geom::Coordinate(prevLon, prevLat));
788 }
789 p_pts->add(geos::geom::Coordinate(lon, lat));
790 prevLon = lon;
791 prevLat = lat;
792 }
793
794
795 // Checks for self-intersection and attempts to correct
796 geos::geom::CoordinateSequence *tempPts = new geos::geom::CoordinateSequence();
797
798 // Gets the starting, second, second to last, and last points to check for validity
799 tempPts->add(geos::geom::Coordinate((*p_pts)[0].x, (*p_pts)[0].y));
800 tempPts->add(geos::geom::Coordinate((*p_pts)[1].x, (*p_pts)[1].y));
801 tempPts->add(geos::geom::Coordinate((*p_pts)[p_pts->size()-3].x, (*p_pts)[p_pts->size()-3].y));
802 tempPts->add(geos::geom::Coordinate((*p_pts)[p_pts->size()-2].x, (*p_pts)[p_pts->size()-2].y));
803 tempPts->add(geos::geom::Coordinate((*p_pts)[0].x, (*p_pts)[0].y));
804
805 geos::geom::Polygon *tempPoly = globalFactory->createPolygon
806 (globalFactory->createLinearRing(*tempPts)).release();
807
808 // Remove the last point of the sequence if it produces invalid polygons
809 if (!tempPoly->isValid()) {
810 std::vector<geos::geom::Coordinate> coords;
811 p_pts->toVector(coords);
812 coords.erase(coords.begin() + p_pts->size() - 2);
813 p_pts->setPoints(coords);
814 }
815
816 delete tempPts;
817 tempPts = NULL;
818 // end self-intersection check
819
820 FixPolePoly(crossingPoints);
821 delete crossingPoints;
822 crossingPoints = NULL;
823 }
824
825
833 void ImagePolygon::FixPolePoly(std::vector<geos::geom::Coordinate> *crossingPoints) {
834 // We currently do not support both poles in one image
835 bool hasNorthPole = false;
836 bool hasSouthPole = false;
837
838 if (p_gMap->SetUniversalGround(90, 0)) {
839 double nPoleSample = Null;
840 double nPoleLine = Null;
841
842 if (p_gMap->Projection()) {
843 nPoleSample = p_gMap->Projection()->WorldX();
844 nPoleLine = p_gMap->Projection()->WorldY();
845 }
846 else if (p_gMap->Camera()) {
847 nPoleSample = p_gMap->Camera()->Sample();
848 nPoleLine = p_gMap->Camera()->Line();
849 }
850
851 if (nPoleSample >= 0.5 && nPoleLine >= 0.5 &&
852 nPoleSample <= p_cube->sampleCount() + 0.5 &&
853 nPoleLine <= p_cube->lineCount() + 0.5 &&
854 SetImage(nPoleSample, nPoleLine)) {
855 hasNorthPole = true;
856 }
857 }
858
859 if (p_gMap->SetUniversalGround(-90, 0)) {
860 double sPoleSample = Null;
861 double sPoleLine = Null;
862
863 if (p_gMap->Projection()) {
864 sPoleSample = p_gMap->Projection()->WorldX();
865 sPoleLine = p_gMap->Projection()->WorldY();
866 }
867 else if (p_gMap->Camera()) {
868 sPoleSample = p_gMap->Camera()->Sample();
869 sPoleLine = p_gMap->Camera()->Line();
870 }
871
872 if (sPoleSample >= 0.5 && sPoleLine >= 0.5 &&
873 sPoleSample <= p_cube->sampleCount() + 0.5 &&
874 sPoleLine <= p_cube->lineCount() + 0.5 &&
875 SetImage(sPoleSample, sPoleLine)) {
876 hasSouthPole = true;
877 }
878 }
879
880 if (hasNorthPole && hasSouthPole) {
881 std::string msg = "Unable to create image footprint because image has both poles";
882 throw IException(IException::Programmer, msg, _FILEINFO_);
883 }
884 else if (crossingPoints->size() == 0) {
885 // No crossing points
886 return;
887 }
888
889 if (hasNorthPole) {
891
892 // If the (north) pole is settable but not within proper angles,
893 // then the polygon does not contain the (north) pole when the cube does
894 if (p_gMap->Camera() &&
896 return;
897 }
898 if (p_gMap->Camera() &&
900 return;
901 }
902 }
903 else if (hasSouthPole) {
904 p_gMap->SetUniversalGround(-90, 0);
905
906 // If the (south) pole is settable but not within proper angles,
907 // then the polygon does not contain the (south) pole when the cube does
908 if (p_gMap->Camera() &&
910 return;
911 }
912 if (p_gMap->Camera() &&
914 return;
915 }
916 }
917
918 geos::geom::Coordinate *closestPoint = &crossingPoints->at(0);
919 geos::geom::Coordinate *pole = NULL;
920
921 // Setup the right pole
922 if (hasNorthPole) {
923 pole = new geos::geom::Coordinate(0, 90);
924 }
925 else if (hasSouthPole) {
926 pole = new geos::geom::Coordinate(0, -90);
927 }
928 else if (crossingPoints->size() % 2 == 1) {
929 geos::geom::Coordinate nPole(0, 90);
930 double nDist = DBL_MAX;
931 geos::geom::Coordinate sPole(0, -90);
932 double sDist = DBL_MAX;
933
934 for (unsigned int index = 0; index < p_pts->size(); index ++) {
935 double north = DistanceSquared(&nPole, &((*p_pts)[index]));
936 if (north < nDist) {
937 nDist = north;
938 }
939 double south = DistanceSquared(&sPole, &((*p_pts)[index]));
940 if (south < sDist) {
941 sDist = south;
942 }
943 }
944
945 if (sDist < nDist) {
946 pole = new geos::geom::Coordinate(0, -90);
947 }
948 else {
949 pole = new geos::geom::Coordinate(0, 90);
950 }
951 }
952
953 // Skip the fix if no pole was determined
954 if (pole == NULL) {
955 return;
956 }
957
958 // Find where the pole needs to be split
959 double closestDistance = DBL_MAX;
960 for (unsigned int i = 0; i < crossingPoints->size(); i++) {
961 geos::geom::Coordinate *temp = &crossingPoints->at(i);
962 double tempDistance;
963 // Make sure the Lat is as close to 0 as possible for a correct distance
964 if (temp->x > 180) {
965 double mod = 0.0;
966 while ((temp->x - mod) > 180) mod += 360.0;
967
968 // Create the modified Point, create a pointer to it, and find the distance
969 geos::geom::Coordinate modPointMem = geos::geom::Coordinate(temp->x - mod, temp->y);
970 geos::geom::Coordinate *modPoint = &modPointMem;
971 tempDistance = DistanceSquared(modPoint, pole);
972 }
973 else {
974 tempDistance = DistanceSquared(temp, pole);
975 }
976 if (tempDistance < closestDistance) {
977 closestDistance = tempDistance;
978 closestPoint = temp;
979 }
980 }
981
982 if (closestDistance == DBL_MAX) {
983 std::string msg = "Image contains a pole but did not detect a meridian crossing!";
984 throw IException(IException::Programmer, msg, _FILEINFO_);
985 }
986
987 // split image at the pole
988 geos::geom::CoordinateSequence *new_points = new geos::geom::CoordinateSequence();
989 for (unsigned int i = 0; i < p_pts->size(); i++) {
990 geos::geom::Coordinate temp = p_pts->getAt(i);
991 new_points->add(temp);
992 if (temp.equals(*closestPoint)) {
993 // Setup direction
994 if (i + 1 != p_pts->size()) { // Prevent overflow exception
995 double fromLon, toLon;
996 if ((p_pts->getAt(i + 1).x - closestPoint->x) > 0) {
997 fromLon = 0.0;
998 toLon = 360.0;
999 }
1000 else {
1001 fromLon = 360.0;
1002 toLon = 0.0;
1003 }
1004
1005 geos::algorithm::LineIntersector lineIntersector;
1006 geos::geom::Coordinate crossingPoint;
1007 geos::geom::Coordinate nPole(0.0, 90.0);
1008 geos::geom::Coordinate sPole(0.0, -90.0);
1009 double dist = DBL_MAX;
1010
1011 for (int num = 0; num < 2 && dist > 180.0; num ++) {
1012 nPole = geos::geom::Coordinate(num * 360.0, 90.0);
1013 sPole = geos::geom::Coordinate(num * 360.0, -90.0);
1014
1015 if (temp.x > 0.0 && p_pts->getAt(i + 1).x > 0.0) {
1016 crossingPoint = geos::geom::Coordinate(p_pts->getAt(i + 1).x - 360.0 + (num * 720.0), p_pts->getAt(i + 1).y);
1017 }
1018 else if (temp.x < 0.0 && p_pts->getAt(i + 1).x < 0.0) { // This should never happen
1019 crossingPoint = geos::geom::Coordinate(p_pts->getAt(i + 1).x + 360.0 - (num * 720.0), p_pts->getAt(i + 1).y);
1020 }
1021
1022 dist = std::sqrt(DistanceSquared(&temp, &crossingPoint));
1023 }
1024
1025 lineIntersector.computeIntersection(nPole, sPole, temp, crossingPoint);
1026
1027 if (lineIntersector.hasIntersection()) {
1028 const geos::geom::Coordinate &intersection = lineIntersector.getIntersection(0);
1029
1030 // Calculate the latituded of the points along the meridian
1031 if (pole->y < intersection.y) {
1032 dist = -dist;
1033 }
1034 vector<double> lats;
1035 double maxLat = std::max(intersection.y, pole->y);
1036 double minLat = std::min(intersection.y, pole->y);
1037 for (double lat = intersection.y + dist; lat < maxLat && lat > minLat; lat += dist) {
1038 lats.push_back(lat);
1039 }
1040
1041 // Add the new points
1042 new_points->add(geos::geom::Coordinate(fromLon, intersection.y));
1043 for (int lat = 0; lat < (int)lats.size(); lat ++) {
1044 new_points->add(geos::geom::Coordinate(fromLon, lats[lat]));
1045 }
1046 new_points->add(geos::geom::Coordinate(fromLon, pole->y));
1047 new_points->add(geos::geom::Coordinate(toLon, pole->y));
1048 for (int lat = lats.size() - 1; lat >= 0; lat --) {
1049 new_points->add(geos::geom::Coordinate(toLon, lats[lat]));
1050 }
1051 new_points->add(geos::geom::Coordinate(toLon, intersection.y));
1052 }
1053 else {
1054 std::string msg = "Image contains a pole but could not determine a meridian crossing!";
1055 throw IException(IException::Programmer, msg, _FILEINFO_);
1056 }
1057
1058 }
1059 }
1060 }
1061 delete p_pts;
1062 p_pts = new_points;
1063 delete pole;
1064 }
1065
1066
1079 bool ImagePolygon::SetImage(const double sample, const double line) {
1080 bool found = false;
1081 if (!p_isProjected) {
1082 found = p_gMap->SetImage(sample, line);
1083 if (!found) {
1084 return false;
1085 }
1086 else {
1087 // Check for valid emission and incidence
1088 try {
1089 if (p_gMap->Camera() &&
1091 return false;
1092 }
1093 if (p_gMap->Camera() &&
1095 return false;
1096 }
1097 }
1098 catch(IException &error) {
1099 }
1100
1109 // Make sure we can go back to image coordinates
1110 // This is done because some camera models due to distortion, get
1111 // a lat/lon for samp/line=1:1, but entering that lat/lon does
1112 // not return samp/line =1:1. Ie. moc WA global images
1113 //double lat = p_gMap->UniversalLatitude();
1114 //double lon = p_gMap->UniversalLongitude();
1115 //return p_gMap->SetUniversalGround(lat,lon);
1116
1117 return found;
1118 }
1119 }
1120 else {
1121 // If projected, make sure the pixel DN is valid before worrying about
1122 // geometry.
1123 p_brick->SetBasePosition((int)sample, (int)line, 1);
1124 p_cube->read(*p_brick);
1125 if (Isis::IsNullPixel((*p_brick)[0])) {
1126 return false;
1127 }
1128 else {
1129 return p_gMap->SetImage(sample, line);
1130 }
1131 }
1132 }
1133
1134
1142 bool convertLon = false;
1143 bool negAdjust = false;
1144 bool newCoords = false; // coordinates have been adjusted
1145 geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateSequence();
1146 double lon, lat;
1147 double lonOffset = 0;
1148 double prevLon = p_pts->getAt(0).x;
1149 double prevLat = p_pts->getAt(0).y;
1150
1151 newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
1152 double dist = 0.0;
1153 for (unsigned int i = 1; i < p_pts->getSize(); i++) {
1154 lon = p_pts->getAt(i).x;
1155 lat = p_pts->getAt(i).y;
1156
1157 // check to see if you just crossed the Meridian
1158 if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
1159 newCoords = true;
1160 // if you were already converting then stop (crossed Meridian even number of times)
1161 if (convertLon) {
1162 convertLon = false;
1163 lonOffset = 0;
1164 }
1165 else { // Need to start converting again, deside how to adjust coordinates
1166 if ((lon - prevLon) > 0) {
1167 lonOffset = -360.;
1168 negAdjust = true;
1169 }
1170 else if ((lon - prevLon) < 0) {
1171 lonOffset = 360.;
1172 negAdjust = false;
1173 }
1174 convertLon = true;
1175 }
1176
1177
1178 }
1179
1180 // Change to a minimum calculation
1181 if (newCoords && dist == 0.0) {
1182 double longitude = (lon + lonOffset) - prevLon;
1183 double latitude = lat - prevLat;
1184 dist = std::sqrt((longitude * longitude) + (latitude * latitude));
1185 }
1186
1187 // add coord
1188 newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));
1189
1190 // set current to old
1191 prevLon = lon;
1192 prevLat = lat;
1193 }
1194
1195 // Nothing was done so return
1196 if (!newCoords) {
1197 geos::geom::Polygon *newPoly = globalFactory->createPolygon
1198 (globalFactory->createLinearRing(*newLonLatPts)).release();
1200 delete newLonLatPts;
1201 return;
1202 }
1203
1204 // bisect into seperate polygons
1205 try {
1206 geos::geom::Polygon *newPoly = globalFactory->createPolygon
1207 (globalFactory->createLinearRing(*newLonLatPts)).release();
1208
1209 geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateSequence();
1210 geos::geom::CoordinateSequence *pts2 = new geos::geom::CoordinateSequence();
1211
1212 // Depending on direction of compensation bound accordingly
1213 //***************************************************
1214
1215 // please verify correct if you change these values
1216 //***************************************************
1217 if (negAdjust) {
1218 pts->add(geos::geom::Coordinate(0., 90.));
1219 pts->add(geos::geom::Coordinate(-360., 90.));
1220 pts->add(geos::geom::Coordinate(-360., -90.));
1221 pts->add(geos::geom::Coordinate(0., -90.));
1222 for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1223 pts->add(geos::geom::Coordinate(0.0, lat));
1224 }
1225 pts->add(geos::geom::Coordinate(0., 90.));
1226 pts2->add(geos::geom::Coordinate(0., 90.));
1227 pts2->add(geos::geom::Coordinate(360., 90.));
1228 pts2->add(geos::geom::Coordinate(360., -90.));
1229 pts2->add(geos::geom::Coordinate(0., -90.));
1230 for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1231 pts2->add(geos::geom::Coordinate(0.0, lat));
1232 }
1233 pts2->add(geos::geom::Coordinate(0., 90.));
1234 }
1235 else {
1236 pts->add(geos::geom::Coordinate(360., 90.));
1237 pts->add(geos::geom::Coordinate(720., 90.));
1238 pts->add(geos::geom::Coordinate(720., -90.));
1239 pts->add(geos::geom::Coordinate(360., -90.));
1240 for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1241 pts->add(geos::geom::Coordinate(360.0, lat));
1242 }
1243 pts->add(geos::geom::Coordinate(360., 90.));
1244 pts2->add(geos::geom::Coordinate(360., 90.));
1245 pts2->add(geos::geom::Coordinate(0., 90.));
1246 pts2->add(geos::geom::Coordinate(0., -90.));
1247 pts2->add(geos::geom::Coordinate(360., -90.));
1248 for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
1249 pts2->add(geos::geom::Coordinate(360.0, lat));
1250 }
1251 pts2->add(geos::geom::Coordinate(360., 90.));
1252 }
1253
1254 geos::geom::Polygon *boundaryPoly = globalFactory->createPolygon
1255 (globalFactory->createLinearRing(*pts)).release();
1256 geos::geom::Polygon *boundaryPoly2 = globalFactory->createPolygon
1257 (globalFactory->createLinearRing(*pts2)).release();
1258 /*------------------------------------------------------------------------
1259 / Intersecting the original polygon (converted coordinates) with the
1260 / boundary polygons will create the multi polygons with the converted coordinates.
1261 / These will need to be converted back to the original coordinates.
1262 /-----------------------------------------------------------------------*/
1263 geos::geom::Geometry *intersection = PolygonTools::Intersect(newPoly, boundaryPoly);
1264 geos::geom::MultiPolygon *convertPoly = PolygonTools::MakeMultiPolygon(intersection);
1265 delete intersection;
1266
1267 intersection = PolygonTools::Intersect(newPoly, boundaryPoly2);
1268 geos::geom::MultiPolygon *convertPoly2 = PolygonTools::MakeMultiPolygon(intersection);
1269 delete intersection;
1270
1271 /*------------------------------------------------------------------------
1272 / Adjust points created in the negative space or >360 space to be back in
1273 / the 0-360 world. This will always only need to be done on convertPoly.
1274 / Then add geometries to finalpolys.
1275 /-----------------------------------------------------------------------*/
1276 std::vector<const geos::geom::Geometry *> *finalpolys = new std::vector<const geos::geom::Geometry *>;
1277 const geos::geom::Geometry *newGeom = NULL;
1278
1279 for (size_t i = 0; i < convertPoly->getNumGeometries(); i++) {
1280 newGeom = (convertPoly->getGeometryN(i))->clone().release();
1281 geos::geom::CoordinateSequence *pts3 = convertPoly->getGeometryN(i)->getCoordinates().release();
1282 geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateSequence();
1283 // fix the points
1284
1285 for (size_t k = 0; k < pts3->getSize() ; k++) {
1286 double lon = pts3->getAt(k).x;
1287 double lat = pts3->getAt(k).y;
1288 if (negAdjust) {
1289 lon = lon + 360;
1290 }
1291 else {
1292 lon = lon - 360;
1293 }
1294 newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);
1295
1296 }
1297 // Add the points to polys
1298 finalpolys->push_back(globalFactory->createPolygon
1299 (globalFactory->createLinearRing(*newLonLatPts)).release());
1300 }
1301
1302 // This loop is over polygons that will always be in 0-360 space no need to convert
1303 for (unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
1304 newGeom = (convertPoly2->getGeometryN(i))->clone().release();
1305 finalpolys->push_back(newGeom);
1306 }
1307
1308 p_polygons = globalFactory->createMultiPolygon(*finalpolys).release();
1309
1310 delete newLonLatPts;
1311 delete pts;
1312 delete pts2;
1313 }
1314 catch(geos::util::IllegalArgumentException *geosIll) {
1315 std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1316 msg += "geos illegal argument [" + IString(geosIll->what()) + "]";
1317 delete geosIll;
1318 throw IException(IException::Unknown, msg, _FILEINFO_);
1319 }
1320 catch(geos::util::GEOSException *geosExc) {
1321 std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1322 msg += "geos exception [" + IString(geosExc->what()) + "]";
1323 delete geosExc;
1324 throw IException(IException::Unknown, msg, _FILEINFO_);
1325 }
1326 catch(IException &e) {
1327 std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1328 msg += "isis operation exception [" + IString(e.what()) + "]";
1329 throw IException(e, IException::Unknown, msg, _FILEINFO_);
1330 }
1331 catch(std::exception &e) {
1332 std::string msg = "Caught std::exception: ";
1333 msg += e.what();
1334 throw IException(IException::Unknown, msg, _FILEINFO_);
1335 }
1336 catch(...) {
1337 std::string msg = "Unable to create image footprint (Fix360Poly) due to ";
1338 msg += "unknown exception";
1339 throw IException(IException::Unknown, msg, _FILEINFO_);
1340 }
1341 }
1342
1343
1352 // Check to see p_polygons is valid data
1353 if (!p_polygons) {
1354 string msg = "Cannot write a NULL polygon!";
1355 throw IException(IException::Programmer, msg, _FILEINFO_);
1356 }
1357
1358 geos::io::WKTWriter *wkt = new geos::io::WKTWriter();
1359 wkt->setTrim(true);
1360
1361 string polyStr = wkt->write(p_polygons);
1362 delete wkt;
1363
1364 Blob newBlob("Footprint", "Polygon");
1365 newBlob.setData(polyStr.c_str(), polyStr.size());
1366 return newBlob;
1367 }
1368
1369
1378 double ImagePolygon::DistanceSquared(const geos::geom::Coordinate *p1, const geos::geom::Coordinate *p2) {
1379 return ((p2->x - p1->x) * (p2->x - p1->x)) + ((p2->y - p1->y) * (p2->y - p1->y));
1380 }
1381
1382
1389 bool hasFourCorners = true;
1390 hasFourCorners &= SetImage(1, 1);
1391 hasFourCorners &= SetImage(p_cubeSamps, 1);
1392 hasFourCorners &= SetImage(p_cubeSamps, p_cubeLines);
1393 hasFourCorners &= SetImage(1, p_cubeLines);
1394 return !hasFourCorners;
1395 }
1396
1397
1398
1412 geos::geom::Coordinate ImagePolygon::FindBestPoint(geos::geom::Coordinate *currentPoint,
1413 geos::geom::Coordinate newPoint,
1414 geos::geom::Coordinate lastPoint) {
1415 geos::geom::Coordinate result = newPoint;
1416 // Use a binary search to snap to the limb when needed
1417 if (p_sampinc > 1 || p_lineinc > 1) {
1418
1419 // Pull the invalid point back inside the image
1420 double x = lastPoint.x;
1421 double y = lastPoint.y;
1422 if (x < p_cubeStartSamp) {
1423 x = p_cubeStartSamp;
1424 }
1425 else if (x > p_cubeSamps) {
1426 x = p_cubeSamps;
1427 }
1428 if (y < p_cubeStartLine) {
1429 y = p_cubeStartLine;
1430 }
1431 else if (y > p_cubeLines) {
1432 y = p_cubeLines;
1433 }
1434 geos::geom::Coordinate invalid(x, y);
1435 geos::geom::Coordinate valid(result.x, result.y);
1436
1437 // Find the best valid Coordinate
1438 while (!SetImage(invalid.x, invalid.y)) {
1439 int x, y;
1440 if (invalid.x > valid.x) {
1441 x = (int)invalid.x - 1;
1442 }
1443 else if (invalid.x < valid.x) {
1444 x = (int)invalid.x + 1;
1445 }
1446 else {
1447 x = (int)invalid.x;
1448 }
1449 if (invalid.y > valid.y) {
1450 y = (int)invalid.y - 1;
1451 }
1452 else if (invalid.y < valid.y) {
1453 y = (int)invalid.y + 1;
1454 }
1455 else {
1456 y = (int)invalid.y;
1457 }
1458 invalid = geos::geom::Coordinate(x, y);
1459 }
1460
1461 result = FixCornerSkip(currentPoint, invalid);
1462 }
1463
1464 return result;
1465 }
1466
1467
1477 geos::geom::Coordinate ImagePolygon::FixCornerSkip(geos::geom::Coordinate *currentPoint,
1478 geos::geom::Coordinate newPoint) {
1479 geos::geom::Coordinate originalPoint = newPoint;
1480 geos::geom::Coordinate modPoint = newPoint;
1481
1482 // Step Too big
1484 return newPoint;
1485 }
1486
1487 // An upper left corner
1488 else if (currentPoint->x < newPoint.x && currentPoint->y > newPoint.y) {
1489 while (newPoint.x >= currentPoint->x && SetImage(newPoint.x, newPoint.y)) {
1490 modPoint = newPoint;
1491 newPoint.x -= 1;
1492 }
1493 }
1494
1495 // An upper right corner
1496 else if (currentPoint->y < newPoint.y && currentPoint->x < newPoint.x) {
1497 while (newPoint.y >= currentPoint->y && SetImage(newPoint.x, newPoint.y)) {
1498 modPoint = newPoint;
1499 newPoint.y -= 1;
1500 }
1501 }
1502
1503 // An lower right corner
1504 else if (currentPoint->x > newPoint.x && currentPoint->y < newPoint.y) {
1505 while (newPoint.x <= currentPoint->x && SetImage(newPoint.x, newPoint.y)) {
1506 modPoint = newPoint;
1507 newPoint.x += 1;
1508 }
1509 }
1510
1511 // An lower left corner
1512 else if (currentPoint->y > newPoint.y && currentPoint->x > newPoint.x) {
1513 while (newPoint.y <= currentPoint->y && SetImage(newPoint.x, newPoint.y)) {
1514 modPoint = newPoint;
1515 newPoint.y += 1;
1516 }
1517 }
1518
1519 if (currentPoint->x == modPoint.x && currentPoint->y == modPoint.y) {
1520 return originalPoint;
1521 }
1522 return modPoint;
1523 }
1524
1525
1533 void ImagePolygon::FindSubpixel(std::vector<geos::geom::Coordinate> & points) {
1534 if (p_subpixelAccuracy > 0) {
1535
1536 // Fix the polygon with subpixel accuracy
1537 geos::geom::Coordinate old = points.at(0);
1538 bool didStartingPoint = false;
1539 for (unsigned int pt = 1; !didStartingPoint; pt ++) {
1540 if (pt >= points.size() - 1) {
1541 pt = 0;
1542 didStartingPoint = true;
1543 }
1544
1545 // Binary Coordinate Search
1546 double maxStep = std::max(p_sampinc, p_lineinc);
1547 double stepY = (old.x - points.at(pt + 1).x) / maxStep;
1548 double stepX = (points.at(pt + 1).y - old.y) / maxStep;
1549
1550 geos::geom::Coordinate valid = points.at(pt);
1551 geos::geom::Coordinate invalid(valid.x + stepX, valid.y + stepY);
1552
1553 for (int itt = 0; itt < p_subpixelAccuracy; itt ++) {
1554 geos::geom::Coordinate half((valid.x + invalid.x) / 2.0, (valid.y + invalid.y) / 2.0);
1555 if (SetImage(half.x, half.y) && InsideImage(half.x, half.y)) {
1556 valid = half;
1557 }
1558 else {
1559 invalid = half;
1560 }
1561 }
1562
1563 old = points.at(pt);
1564
1565 // Set new coordinate
1566 points[pt] = valid;
1567
1568 }
1569
1570 // Fix starting point
1571 points[points.size()-1] = geos::geom::Coordinate(points[0].x, points[0].y);
1572
1573 }
1574 }
1575
1576
1577} // end namespace isis
Buffer for containing a three dimensional section of an image.
Definition Brick.h:45
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:120
virtual double Line() const
Returns the current line number.
Definition Camera.cpp:2740
virtual double Sample() const
Returns the current sample number.
Definition Camera.cpp:2720
IO Handler for Isis Cubes.
Definition Cube.h:168
int lineCount() const
Definition Cube.cpp:1740
Camera * camera()
Return a camera associated with the cube.
Definition Cube.cpp:1457
int sampleCount() const
Definition Cube.cpp:1813
PixelType pixelType() const
Definition Cube.cpp:1764
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:813
virtual QString fileName() const
Returns the opened cube's filename.
Definition Cube.cpp:1569
Projection * projection()
Definition Cube.cpp:1800
Isis exception class.
Definition IException.h:91
@ Unknown
A type of error that cannot be classified as any of the other error types.
Definition IException.h:118
@ User
A type of error that could only have occurred due to a mistake on the user's part (e....
Definition IException.h:126
@ Programmer
This error is for when a programmer made an API call that was illegal.
Definition IException.h:146
Adds specific functionality to C++ strings.
Definition IString.h:165
Create cube polygons, read/write polygons to blobs.
bool IsLimb()
Returns True when the input image is a limb image.
double validSampleDim()
Retuns the maximum valid sample width of the cube set with either initCube() or Create()
geos::geom::Coordinate FindFirstPoint()
Finds the first point that projects in an image.
Brick * p_brick
Used to check for valid DNs.
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.
Cube * p_cube
The cube provided.
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.
int p_sampinc
The increment for walking along the polygon in the sample direction.
ImagePolygon()
Constructs a Polygon object, setting the polygon name.
int p_cubeLines
The number of lines in the cube.
geos::geom::CoordinateSequence * p_pts
The sequence of coordinates that compose the boundary of the image.
int p_cubeSamps
The number of samples in the cube.
int p_subpixelAccuracy
The subpixel accuracy to use.
int p_lineinc
The increment for walking along the polygon in the line direcction.
geos::geom::Coordinate * m_botCoord
The cube's bot-most valid coord.
geos::geom::Coordinate * m_leftCoord
The cube's left-most valid coord.
UniversalGroundMap * p_gMap
The cube's ground map.
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...
int p_cubeStartLine
The the line of the first valid point in the cube.
double p_incidence
The maximum incidence angle to consider valid.
double DistanceSquared(const geos::geom::Coordinate *p1, const geos::geom::Coordinate *p2)
Calculates the distance squared between two coordinates.
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.
void FindSubpixel(std::vector< geos::geom::Coordinate > &points)
Takes p_polygons in sample/line space and finds its subpixel accuracy.
std::string p_polyStr
The string representation of the polygon.
bool p_ellipsoid
Uses an ellipsoid if a limb is detected.
bool SetImage(const double sample, const double line)
Sets the sample/line values of the cube to get lat/lon values.
void Fix360Poly()
If the cube crosses the 0/360 boundary and does not include a pole, the polygon is separated into mul...
~ImagePolygon()
Destroys the Polygon object.
void WalkPoly()
Walks the image finding its lon lat polygon and stores it to p_pts.
std::string polyStr() const
Return a geos Multipolygon.
double validLineDim()
Retuns the maximum valid line width of the cube set with either initCube() or Create()
bool p_isProjected
True when the provided cube is projected.
bool InsideImage(double sample, double line)
This returns true if sample/line are inside the cube.
void calcImageBorderCoordinates()
Calculates the four border points, particularly for validSampleDim() and validLineDim()
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.
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...
geos::geom::Coordinate * m_topCoord
The cube's top-most valid coord.
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...
Blob toBlob() const
Serialize the ImagePolygon to a Blob.
geos::geom::Coordinate * m_rightCoord
The cube's right-most valid coord.
geos::geom::MultiPolygon * p_polygons
The multipolygon of the image.
int p_cubeStartSamp
The the sample of the first valid point in the cube.
static geos::geom::MultiPolygon * MakeMultiPolygon(const geos::geom::Geometry *geom)
Make a geos::geom::MultiPolygon out of the components of the argument.
static geos::geom::Geometry * Intersect(const geos::geom::Geometry *geom1, const geos::geom::Geometry *geom2)
This applies the geos Intersect operator.
virtual double WorldY() const
This returns the world Y coordinate provided SetGround, SetCoordinate, SetUniversalGround,...
virtual double WorldX() const
This returns the world X coordinate provided SetGround, SetCoordinate, SetUniversalGround,...
Container for cube-like labels.
Definition Pvl.h:119
void IgnoreElevationModel(bool ignore)
This allows you to ignore the cube elevation model and use the ellipse.
Definition Sensor.cpp:62
virtual double IncidenceAngle() const
Returns the incidence angle in degrees.
Definition Sensor.cpp:339
virtual double EmissionAngle() const
Returns the emission angle in degrees.
Definition Sensor.cpp:328
Universal Ground Map.
double UniversalLongitude() const
Returns the universal longitude of the camera model or projection.
bool SetImage(double sample, double line)
Returns whether the sample/line postion was set successfully in the camera model or projection.
void SetBand(const int band)
Set the image band number.
double UniversalLatitude() const
Returns the universal latitude of the camera model or projection.
bool SetUniversalGround(double lat, double lon)
Returns whether the lat/lon position was set successfully in the camera model or projection.
Isis::Projection * Projection() const
Return the projection associated with the ground map (NULL implies none)
Isis::Camera * Camera() const
Return the camera associated with the ground map (NULL implies none)
This is free and unencumbered software released into the public domain.
Definition Apollo.h:16
bool IsNullPixel(const double d)
Returns if the input pixel is null.
const double Null
Value for an Isis Null pixel.
Namespace for the standard library.