24 #include "IsisDebug.h"
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>
56 ImagePolygon::ImagePolygon() :
Blob(
"Footprint",
"Polygon") {
115 int ns,
int nl,
int band) {
133 QString msg =
"Can not create polygon, ";
135 msg +=
"] is not a camera or map projection";
139 polyError.
append(camError);
140 polyError.
append(projError);
157 p_cubeSamps = std::min(p_cubeSamps, ss + ns);
172 std::string msg =
"Cannot use an ellipsoid shape model";
173 msg +=
" on a limb image without a camera.";
211 int ss,
int sl,
int ns,
int nl,
int band,
212 bool increasePrecision) {
216 if (sinc < 1 || linc < 1) {
217 std::string msg =
"Sample and line increments must be 1 or greater";
221 cam =
initCube(cube, ss, sl, ns, nl, band);
224 bool polygonGenerated =
false;
225 while (!polygonGenerated) {
231 p_pts =
new geos::geom::CoordinateArraySequence();
235 polygonGenerated =
true;
243 if (increasePrecision && (sinc > 1 || linc > 1)) {
247 QString msg =
"Cannot find polygon for image "
249 msg += increasePrecision ?
"Cannot increase precision any further" :
250 "The increment/step size might be too large";
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;
299 if (recursionDepth > 6) {
300 return *currentPoint;
304 if (x == 0.0 && y == 0.0) {
307 double s = currentPoint->x + samp;
308 double l = currentPoint->y + line;
312 geos::geom::Coordinate next(s, l);
318 std::string msg =
"Unable to create image footprint. Starting point is not on the edge of the image.";
321 else if (x < 0 && y < 0) {
322 geos::geom::Coordinate next(currentPoint->x, currentPoint->y - 1 *
p_lineinc);
325 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
331 else if (x == 0.0 && y < 0) {
332 geos::geom::Coordinate next(currentPoint->x + 1 *
p_sampinc, currentPoint->y - 1 *
p_lineinc);
335 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
341 else if (x > 0 && y < 0) {
342 geos::geom::Coordinate next(currentPoint->x + 1 *
p_sampinc, currentPoint->y);
345 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
351 else if (x > 0 && y == 0.0) {
352 geos::geom::Coordinate next(currentPoint->x + 1 *
p_sampinc, currentPoint->y + 1 *
p_lineinc);
355 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
361 else if (x > 0 && y > 0) {
362 geos::geom::Coordinate next(currentPoint->x, currentPoint->y + 1 *
p_lineinc);
365 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
371 else if (x == 0.0 && y > 0) {
372 geos::geom::Coordinate next(currentPoint->x - 1 *
p_sampinc, currentPoint->y + 1 *
p_lineinc);
375 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
381 else if (x < 0 && y > 0) {
382 geos::geom::Coordinate next(currentPoint->x - 1 *
p_sampinc, currentPoint->y);
385 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
391 else if (x < 0 && y == 0.0) {
392 geos::geom::Coordinate next(currentPoint->x - 1 *
p_sampinc, currentPoint->y - 1 *
p_lineinc);
395 result =
FindNextPoint(currentPoint, next, recursionDepth + 1);
402 std::string msg =
"Unable to create image footprint. Error walking image.";
432 const double origSample = sample - sinc;
434 const double origLine = line - linc;
437 if (sample < startSample && sinc < 0) {
439 if (origSample == startSample) {
443 sample = startSample;
447 if (sample > endSample && sinc > 0) {
449 if (origSample == endSample) {
457 if (line < startLine && linc < 0) {
459 if (origLine == startLine) {
467 if (line > endLine && linc > 0) {
469 if (fabs(origLine - endLine) < 0.5) {
510 geos::geom::Coordinate firstPoint(sample, line);
511 geos::geom::Coordinate lastPoint = firstPoint;
512 if (!firstPoint.equals(
FindNextPoint(&firstPoint, lastPoint))) {
520 std::string msg =
"No lat/lon data found for image";
569 m_leftCoord =
new geos::geom::Coordinate(sample, line);
578 m_rightCoord =
new geos::geom::Coordinate(sample, line);
588 m_topCoord =
new geos::geom::Coordinate(sample, line);
598 m_botCoord =
new geos::geom::Coordinate(sample, line);
614 vector<geos::geom::Coordinate> points;
615 double lat, lon, prevLat, prevLon;
619 points.push_back(firstPoint);
625 geos::geom::Coordinate currentPoint = firstPoint;
626 geos::geom::Coordinate lastPoint = firstPoint;
627 geos::geom::Coordinate tempPoint;
635 bool snapToFirstPoint =
true;
641 snapToFirstPoint &= (points.size() > 2);
648 snapToFirstPoint &= (
DistanceSquared(¤tPoint, &firstPoint) < (minStepSize * minStepSize));
651 if (snapToFirstPoint) {
652 tempPoint = firstPoint;
658 for (
int pt = 0; pt < (int)points.size(); pt ++) {
659 if (points[pt].equals(tempPoint)) {
660 tempPoint = firstPoint;
667 if (tempPoint.equals(currentPoint)) {
668 geos::geom::Coordinate oldDuplicatePoint = tempPoint;
671 tempPoint = lastPoint;
672 lastPoint = currentPoint;
673 currentPoint = tempPoint;
677 if (points.size() < 3) {
678 std::string msg =
"Failed to find next point in the image.";
687 if (tempPoint.equals(currentPoint) || tempPoint.equals(oldDuplicatePoint)) {
688 std::string msg =
"Failed to find next valid point in the image.";
696 if (points[points.size()-3].x == tempPoint.x &&
697 points[points.size()-3].y == tempPoint.y) {
703 currentPoint = points[points.size()-1];
715 if (points.size() > 250) {
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]) {
729 if (cycleStart != 0) {
730 vector<geos::geom::Coordinate> cyclePoints;
731 for (
int pt = cycleStart; pt <= cycleEnd; pt ++) {
732 cyclePoints.push_back(points[pt]);
735 points = cyclePoints;
742 lastPoint = currentPoint;
743 currentPoint = tempPoint;
744 points.push_back(currentPoint);
747 while (!currentPoint.equals(firstPoint));
749 if (points.size() <= 3) {
750 std::string msg =
"Failed to find enough points on the image.";
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));
766 if (abs(lon - prevLon) >= 180 && i != 0) {
767 crossingPoints->push_back(geos::geom::Coordinate(prevLon, prevLat));
769 p_pts->add(geos::geom::Coordinate(lon, lat));
776 geos::geom::CoordinateSequence *tempPts =
new geos::geom::CoordinateArraySequence();
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));
785 geos::geom::Polygon *tempPoly = globalFactory.createPolygon
786 (globalFactory.createLinearRing(tempPts), NULL);
789 if (!tempPoly->isValid()) {
798 delete crossingPoints;
799 crossingPoints = NULL;
812 bool hasNorthPole =
false;
813 bool hasSouthPole =
false;
816 double nPoleSample =
Null;
817 double nPoleLine =
Null;
828 if (nPoleSample >= 0.5 && nPoleLine >= 0.5 &&
829 nPoleSample <= p_cube->sampleCount() + 0.5 &&
830 nPoleLine <= p_cube->lineCount() + 0.5) {
836 double sPoleSample =
Null;
837 double sPoleLine =
Null;
848 if (sPoleSample >= 0.5 && sPoleLine >= 0.5 &&
849 sPoleSample <= p_cube->sampleCount() + 0.5 &&
850 sPoleLine <= p_cube->lineCount() + 0.5) {
855 if (hasNorthPole && hasSouthPole) {
856 std::string msg =
"Unable to create image footprint because image has both poles";
859 else if (crossingPoints->size() == 0) {
878 else if (hasSouthPole) {
893 geos::geom::Coordinate *closestPoint = &crossingPoints->at(0);
894 geos::geom::Coordinate *pole = NULL;
898 pole =
new geos::geom::Coordinate(0, 90);
900 else if (hasSouthPole) {
901 pole =
new geos::geom::Coordinate(0, -90);
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;
909 for (
unsigned int index = 0; index <
p_pts->size(); index ++) {
921 pole =
new geos::geom::Coordinate(0, -90);
924 pole =
new geos::geom::Coordinate(0, 90);
934 double closestDistance = DBL_MAX;
935 for (
unsigned int i = 0; i < crossingPoints->size(); i++) {
936 geos::geom::Coordinate *temp = &crossingPoints->at(i);
941 while ((temp->x - mod) > 180) mod += 360.0;
944 geos::geom::Coordinate modPointMem = geos::geom::Coordinate(temp->x - mod, temp->y);
945 geos::geom::Coordinate *modPoint = &modPointMem;
951 if (tempDistance < closestDistance) {
952 closestDistance = tempDistance;
957 if (closestDistance == DBL_MAX) {
958 std::string msg =
"Image contains a pole but did not detect a meridian crossing!";
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)) {
969 if (i + 1 !=
p_pts->size()) {
970 double fromLon, toLon;
971 if ((
p_pts->getAt(i + 1).x - closestPoint->x) > 0) {
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;
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);
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);
993 else if (temp.x < 0.0 &&
p_pts->getAt(i + 1).x < 0.0) {
994 crossingPoint = geos::geom::Coordinate(
p_pts->getAt(i + 1).x + 360.0 - (num * 720.0),
p_pts->getAt(i + 1).y);
1000 lineIntersector.computeIntersection(nPole, sPole, temp, crossingPoint);
1002 if (lineIntersector.hasIntersection()) {
1003 const geos::geom::Coordinate &intersection = lineIntersector.getIntersection(0);
1006 if (pole->y < intersection.y) {
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);
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]));
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]));
1026 new_points->add(geos::geom::Coordinate(toLon, intersection.y));
1029 std::string msg =
"Image contains a pole but could not determine a meridian crossing!";
1117 bool convertLon =
false;
1118 bool negAdjust =
false;
1119 bool newCoords =
false;
1120 geos::geom::CoordinateSequence *newLonLatPts =
new geos::geom::CoordinateArraySequence();
1122 double lonOffset = 0;
1123 double prevLon =
p_pts->getAt(0).x;
1124 double prevLat =
p_pts->getAt(0).y;
1126 newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
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;
1132 if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
1140 if ((lon - prevLon) > 0) {
1144 else if ((lon - prevLon) < 0) {
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));
1162 newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));
1171 geos::geom::Polygon *newPoly = globalFactory.createPolygon
1172 (globalFactory.createLinearRing(newLonLatPts), NULL);
1174 delete newLonLatPts;
1180 geos::geom::Polygon *newPoly = globalFactory.createPolygon
1181 (globalFactory.createLinearRing(newLonLatPts), NULL);
1183 geos::geom::CoordinateSequence *pts =
new geos::geom::CoordinateArraySequence();
1184 geos::geom::CoordinateSequence *pts2 =
new geos::geom::CoordinateArraySequence();
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));
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));
1207 pts2->add(geos::geom::Coordinate(0., 90.));
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));
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));
1225 pts2->add(geos::geom::Coordinate(360., 90.));
1228 geos::geom::Polygon *boundaryPoly = globalFactory.createPolygon
1229 (globalFactory.createLinearRing(pts), NULL);
1230 geos::geom::Polygon *boundaryPoly2 = globalFactory.createPolygon
1231 (globalFactory.createLinearRing(pts2), NULL);
1239 delete intersection;
1243 delete intersection;
1250 vector<geos::geom::Geometry *> *finalpolys =
new vector<geos::geom::Geometry *>;
1251 geos::geom::Geometry *newGeom = NULL;
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();
1259 for (
unsigned int k = 0; k < pts->getSize() ; k++) {
1260 double lon = pts->getAt(k).x;
1261 double lat = pts->getAt(k).y;
1268 newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);
1272 finalpolys->push_back(globalFactory.createPolygon
1273 (globalFactory.createLinearRing(newLonLatPts), NULL));
1277 for (
unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
1278 newGeom = (convertPoly2->getGeometryN(i))->clone();
1279 finalpolys->push_back(newGeom);
1282 p_polygons = globalFactory.createMultiPolygon(*finalpolys);
1286 delete newLonLatPts;
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()) +
"]";
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()) +
"]";
1303 std::string msg =
"Unable to create image footprint (Fix360Poly) due to ";
1304 msg +=
"isis operation exception [" +
IString(e.
what()) +
"]";
1308 std::string msg =
"Unable to create image footprint (Fix360Poly) due to ";
1309 msg +=
"unknown exception";
1326 is.seekg(sbyte, std::ios::beg);
1328 QString msg =
"Error preparing to read data from " +
p_type +
1343 QString msg =
"Error reading data from " +
p_type +
" [" +
1348 geos::io::WKTReader *wkt =
new geos::io::WKTReader(&(globalFactory));
1356 geos::io::WKTWriter *wkt =
new geos::io::WKTWriter();
1360 string msg =
"Cannot write a NULL polygon!";
1391 return ((p2->x - p1->x) * (p2->x - p1->x)) + ((p2->y - p1->y) * (p2->y - p1->y));
1401 bool hasFourCorners =
true;
1406 return !hasFourCorners;
1425 geos::geom::Coordinate newPoint,
1426 geos::geom::Coordinate lastPoint) {
1427 geos::geom::Coordinate result = newPoint;
1432 double x = lastPoint.x;
1433 double y = lastPoint.y;
1446 geos::geom::Coordinate invalid(x, y);
1447 geos::geom::Coordinate valid(result.x, result.y);
1450 while (!
SetImage(invalid.x, invalid.y)) {
1452 if (invalid.x > valid.x) {
1453 x = (int)invalid.x - 1;
1455 else if (invalid.x < valid.x) {
1456 x = (int)invalid.x + 1;
1461 if (invalid.y > valid.y) {
1462 y = (int)invalid.y - 1;
1464 else if (invalid.y < valid.y) {
1465 y = (int)invalid.y + 1;
1470 invalid = geos::geom::Coordinate(x, y);
1490 geos::geom::Coordinate newPoint) {
1491 geos::geom::Coordinate originalPoint = newPoint;
1492 geos::geom::Coordinate modPoint = newPoint;
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;
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;
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;
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;
1531 if (currentPoint->x == modPoint.x && currentPoint->y == modPoint.y) {
1532 return originalPoint;
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) {
1554 didStartingPoint =
true;
1559 double stepY = (old.x - points.at(pt + 1).x) / maxStep;
1560 double stepX = (points.at(pt + 1).y - old.y) / maxStep;
1562 geos::geom::Coordinate valid = points.at(pt);
1563 geos::geom::Coordinate invalid(valid.x + stepX, valid.y + stepY);
1566 geos::geom::Coordinate half((valid.x + invalid.x) / 2.0, (valid.y + invalid.y) / 2.0);
1575 old = points.at(pt);
1583 points[points.size()-1] = geos::geom::Coordinate(points[0].x, points[0].y);