PolygonTools.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include <string>
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <cmath>
 
#include <QDebug>
 
#include "geos/geom/CoordinateSequence.h"
#include "geos/geom/LinearRing.h"
#include "geos/geom/Point.h"
#include "geos/geom/Polygon.h"
#include "geos/operation/distance/DistanceOp.h"
#include "geos/operation/overlay/snap/SnapOverlayOp.h"
#include "geos/operation/overlay/snap/GeometrySnapper.h"
 
#include "SpecialPixel.h"
#include "PolygonTools.h"
#include "TProjection.h"
#include "ProjectionFactory.h"
#include "UniversalGroundMap.h"
 
using namespace std;
namespace Isis {
 
  geos::geom::MultiPolygon *PolygonTools::LatLonToXY(
      const geos::geom::MultiPolygon &lonLatPolygon, TProjection *projection) {
    if (projection == NULL) {
      string msg = "Unable to convert Lon/Lat polygon to X/Y. ";
      msg += "No projection has was supplied";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    // Convert the Lat/Lon poly coordinates to X/Y coordinates
    if(lonLatPolygon.isEmpty()) {
      return globalFactory->createMultiPolygon().release();
    }
    else {
      vector<const geos::geom::Geometry *> *xyPolys = new vector<const geos::geom::Geometry *>;
      // Convert each polygon in this multi-polygon
      for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); ++g) {
        const geos::geom::Polygon *poly =
            dynamic_cast<const geos::geom::Polygon *>(
              lonLatPolygon.getGeometryN(g));
 
        // Convert each hole inside this polygon
        std::vector<std::unique_ptr<geos::geom::LinearRing>> holes;
        for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
          geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateSequence();
          geos::geom::CoordinateSequence *llcoords =
            poly->getInteriorRingN(h)->getCoordinates().release();
 
          // Convert each coordinate in this hole
          for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
            projection->SetGround(llcoords->getAt(cord).y,
                                  llcoords->getAt(cord).x);
            xycoords->add(geos::geom::Coordinate(projection->XCoord(),
                                                 projection->YCoord()));
          } // end num coords in hole loop
 
          std::unique_ptr<geos::geom::LinearRing> hole = globalFactory->createLinearRing(*xycoords);
 
          if(hole->isValid() && !hole->isEmpty()) {
            holes.push_back(hole->clone());
          }
          else {
            hole.reset();
          }
        } // end num holes in polygon loop
 
        // Convert the exterior ring of this polygon
        geos::geom::CoordinateSequence *xycoords = new geos::geom::CoordinateSequence();
        geos::geom::CoordinateSequence *llcoords =
          poly->getExteriorRing()->getCoordinates().release();
 
        // Convert each coordinate in the exterior ring of this polygon
        for (unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
          projection->SetGround(llcoords->getAt(cord).y,
                                llcoords->getAt(cord).x);
          xycoords->add(geos::geom::Coordinate(projection->XCoord(),
                                               projection->YCoord()));
        } // end exterior ring coordinate loop
 
        geos::geom::Polygon *newPoly = globalFactory->createPolygon(
                                         globalFactory->createLinearRing(*xycoords), std::move(holes)).release();
 
        if(newPoly->isValid() && !newPoly->isEmpty() && newPoly->getArea() > 1.0e-14) {
          xyPolys->push_back(newPoly);
        }
        else {
          delete newPoly;
        }
      } // end num geometry in multi-poly
 
      // Create a new multipoly from all the new X/Y polygon(s)
      geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(*xyPolys).release();
 
      if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
        return spikedPoly;
      }
      else {
        try {
          geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
 
          delete spikedPoly;
          spikedPoly = NULL;
 
          return despikedPoly;
        }
        catch(IException &e) {
          IString msg = "Unable to convert polygon from Lat/Lon to X/Y";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }
 
    } // end else
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::XYToLatLon(
    const geos::geom::MultiPolygon &xYPolygon, TProjection *projection) {
 
    if(projection == NULL) {
      string msg = "Unable to convert X/Y polygon to Lon/Lat. ";
      msg += "No projection was supplied";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    // Convert the X/Y poly coordinates to Lat/Lon coordinates
    if(xYPolygon.isEmpty()) {
      return globalFactory->createMultiPolygon().release();
    }
    else {
      vector<const geos::geom::Geometry *> *llPolys = new vector<const geos::geom::Geometry *>;
      // Convert each polygon in this multi-polygon
      for(unsigned int g = 0; g < xYPolygon.getNumGeometries(); ++g) {
        const geos::geom::Polygon *poly =
            dynamic_cast<const geos::geom::Polygon *>(
              xYPolygon.getGeometryN(g));
 
        // Convert each hole inside this polygon
        std::vector<std::unique_ptr<geos::geom::LinearRing>> holes;
        for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
          geos::geom::CoordinateSequence *llcoords = new geos::geom::CoordinateSequence();
          geos::geom::CoordinateSequence *xycoords =
            poly->getInteriorRingN(h)->getCoordinates().release();
 
          // Convert each coordinate in this hole
          for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
            projection->SetWorld(xycoords->getAt(cord).x,
                                 xycoords->getAt(cord).y);
            llcoords->add(geos::geom::Coordinate(projection->Longitude(),
                                                 projection->Latitude()));
          } // end num coords in hole loop
          holes.push_back(globalFactory->createLinearRing(*llcoords));
        } // end num holes in polygon loop
 
        // Convert the exterior ring of this polygon
        geos::geom::CoordinateSequence *llcoords = new geos::geom::CoordinateSequence();
        geos::geom::CoordinateSequence *xycoords =
          poly->getExteriorRing()->getCoordinates().release();
 
        // Convert each coordinate in the exterior ring of this polygon
        for(unsigned int cord = 0; cord < xycoords->getSize(); ++cord) {
          projection->SetWorld(xycoords->getAt(cord).x,
                               xycoords->getAt(cord).y);
          llcoords->add(geos::geom::Coordinate(projection->Longitude(),
                                               projection->Latitude()));
        } // end exterior ring coordinate loop
 
        llPolys->push_back(globalFactory->createPolygon(
                             globalFactory->createLinearRing(*llcoords), std::move(holes)).release());
      } // end num geometry in multi-poly
 
 
      // Create a new multipoly from all the new Lat/Lon polygon(s)
      geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(*llPolys).release();
 
      if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
        return spikedPoly;
      }
      else {
        try {
          geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
 
          delete spikedPoly;
          spikedPoly = NULL;
 
          return despikedPoly;
        }
        catch(IException &e) {
          IString msg = "Unable to convert polygon from X/Y to Lat/Lon";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }
    } // end else
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::LatLonToSampleLine(
    const geos::geom::MultiPolygon &lonLatPolygon, UniversalGroundMap *ugm) {
 
    if(ugm == NULL) {
      string msg = "Unable to convert Lon/Lat polygon to Sample/Line. ";
      msg += "No UniversalGroundMap was supplied";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    // Convert the Lon/Lat poly coordinates to Sample/Line coordinates
    if(lonLatPolygon.isEmpty()) {
      return globalFactory->createMultiPolygon().release();
    }
    else {
      vector<const geos::geom::Geometry *> slPolys;
      // Convert each polygon in this multi-polygon
      for(unsigned int g = 0; g < lonLatPolygon.getNumGeometries(); g++) {
        const geos::geom::Polygon *poly =
            dynamic_cast<const geos::geom::Polygon *>(
              lonLatPolygon.getGeometryN(g));
 
        // Convert each hole inside this polygon
        vector<unique_ptr<geos::geom::LinearRing>> holes;
        for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
          geos::geom::CoordinateSequence *slcoords = new geos::geom::CoordinateSequence();
          geos::geom::CoordinateSequence *llcoords =
            poly->getInteriorRingN(h)->getCoordinates().release();
 
          // Convert each coordinate in this hole
          for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
            ugm->SetUniversalGround(llcoords->getAt(cord).y,
                                    llcoords->getAt(cord).x);
            slcoords->add(geos::geom::Coordinate(ugm->Sample(),
                                                 ugm->Line()));
          } // end num coords in hole loop
          holes.push_back(globalFactory->createLinearRing(*slcoords));
          delete slcoords;
          delete llcoords;
        } // end num holes in polygon loop
 
        // Convert the exterior ring of this polygon
        geos::geom::CoordinateSequence *slcoords = new geos::geom::CoordinateSequence();
        geos::geom::CoordinateSequence *llcoords =
          poly->getExteriorRing()->getCoordinates().release();
 
        // Convert each coordinate in the exterior ring of this polygon
        for(unsigned int cord = 0; cord < llcoords->getSize(); ++cord) {
          if (ugm->SetUniversalGround(llcoords->getAt(cord).y,
                                      llcoords->getAt(cord).x)) {
            slcoords->add(geos::geom::Coordinate(ugm->Sample(),
                                                ugm->Line()));
          }
        } // end exterior ring coordinate loop
 
        // Make sure that the line string is closed.
        if (slcoords->getSize() > 0 && !slcoords->front().equals(slcoords->back())) {
          slcoords->add(slcoords->front());
        }
 
        try {
          slPolys.push_back(globalFactory->createPolygon(
                            globalFactory->createLinearRing(*slcoords), std::move(holes)).release());
        }
        catch (std::exception &e) {
          throw IException(IException::Unknown,
                           QObject::tr("Unable to convert polygon from Lat/Lon to Sample/Line. The "
                                       "error given was [%1].").arg(e.what()),
                           _FILEINFO_);
        }
 
        delete llcoords;
      } // end num geometry in multi-poly
 
      // Create a new multipoly from all the new Sample/Line polygon(s)
      geos::geom::MultiPolygon *spikedPoly = globalFactory->createMultiPolygon(slPolys).release();
 
      if(spikedPoly->isValid() && !spikedPoly->isEmpty()) {
        return spikedPoly;
      }
      else {
        try {
          geos::geom::MultiPolygon *despikedPoly = Despike(spikedPoly);
 
          delete spikedPoly;
          spikedPoly = NULL;
 
          return despikedPoly;
        }
        catch (IException &e) {
          IString msg = "Unable to convert polygon from Lat/Lon to Sample/Line";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }
    } // end else
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon *mpolygon) {
 
    vector<const geos::geom::Geometry *> *polys = new vector<const geos::geom::Geometry *>;
    for(unsigned int i = 0; i < mpolygon->getNumGeometries(); ++i) {
      polys->push_back((mpolygon->getGeometryN(i))->clone().release());
    }
    return globalFactory->createMultiPolygon(*polys).release();
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::CopyMultiPolygon(const geos::geom::MultiPolygon &mpolygon) {
 
    vector<const geos::geom::Geometry *> *polys = new vector<const geos::geom::Geometry *>;
    for(unsigned int i = 0; i < mpolygon.getNumGeometries(); ++i) {
      polys->push_back((mpolygon.getGeometryN(i))->clone().release());
    }
    return globalFactory->createMultiPolygon(*polys).release();
  }
 
 
  QString PolygonTools::ToGML(const geos::geom::MultiPolygon *mpolygon, QString idString,
                              QString schema) {
 
    ostringstream os;
 
    //Write out the GML header
    os << "<?xml version=\"1.0\" encoding=\"utf-8\" ?>" << endl;
    os << "<ogr:FeatureCollection" << endl;
    os << "    xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"" << endl;
    os << "    xsi:schemaLocation=\"http://ogr.maptools.org/ " << schema << "\"" << endl;
    os << "    xmlns:ogr=\"http://ogr.maptools.org/\"" << endl;
    os << "    xmlns:gml=\"http://www.opengis.net/gml\">" << endl;
    os << "  <gml:boundedBy>" << endl;
    os << "    <gml:Box>" << endl;
    os << "      <gml:coord>";
    os <<          "<gml:X>" <<
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMinX() <<
                       "</gml:X>";
    os <<          "<gml:Y>" <<
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMinY() <<
                       "</gml:Y>";
    os <<      "</gml:coord>" << endl;
    os << "      <gml:coord>";
    os <<          "<gml:X>" <<
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxX() <<
                       "</gml:X>";
    os <<          "<gml:Y>" <<
                       setprecision(15) << mpolygon->getEnvelopeInternal()->getMaxY() <<
                       "</gml:Y>";
    os <<        "</gml:coord>" << endl;
    os << "    </gml:Box>" << endl;
    os << "  </gml:boundedBy>" << endl << endl;
    os << "  <gml:featureMember>" << endl;
    os << "   <ogr:multi_polygon fid=\"0\">" << endl;
    os << "      <ogr:ID>" << idString << "</ogr:ID>" << endl;
    os << "      <ogr:geometryProperty><gml:MultiPolygon><gml:polygonMember>" <<
                     "<gml:Polygon><gml:outerBoundaryIs><gml:LinearRing><gml:coordinates>";
 
 
    for(unsigned int polyN = 0; polyN < mpolygon->getNumGeometries(); polyN++) {
      geos::geom::CoordinateSequence *pts = mpolygon->getGeometryN(polyN)->getCoordinates().release();
 
      for(unsigned int i = 0; i < pts->getSize(); i++) {
        double lon = pts->getAt(i).x;
        double lat = pts->getAt(i).y;
 
        os << setprecision(15) << lon << "," << setprecision(15) << lat << " ";
      }
    }
 
    os << "</gml:coordinates></gml:LinearRing></gml:outerBoundaryIs>" <<
              "</gml:Polygon></gml:polygonMember></gml:MultiPolygon>" <<
              "</ogr:geometryProperty>" << endl;
    os << "    </ogr:multi_polygon>" << endl;
    os << "  </gml:featureMember>" << endl;
    os << "</ogr:FeatureCollection>";
 
    return os.str().c_str();
  }
 
 
  QString PolygonTools::GMLSchema() {
 
    ostringstream os;
 
    os << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl;
    os << "<xs:schema targetNamespace=\"http://ogr.maptools.org/\" "
              "xmlns:ogr=\"http://ogr.maptools.org/\" "
              "xmlns:xs=\"http://www.w3.org/2001/XMLSchema\" "
              "xmlns:gml=\"http://www.opengis.net/gml\" "
              "elementFormDefault=\"qualified\" "
              "version=\"1.0\">" << endl;
    os << "  <xs:import namespace=\"http://www.opengis.net/gml\" "
              "schemaLocation=\"http://schemas.opengis.net/gml/2.1.2/feature.xsd\"/>" << endl;
    os << "  <xs:element name=\"FeatureCollection\" "
                 "type=\"ogr:FeatureCollectionType\" "
                 "substitutionGroup=\"gml:_FeatureCollection\"/>" << endl;
    os << "  <xs:complexType name=\"FeatureCollectionType\">" << endl;
    os << "    <xs:complexContent>" << endl;
    os << "      <xs:extension base=\"gml:AbstractFeatureCollectionType\">" << endl;
    os << "        <xs:attribute name=\"lockId\" type=\"xs:string\" use=\"optional\"/>" << endl;
    os << "        <xs:attribute name=\"scope\" type=\"xs:string\" use=\"optional\"/>" << endl;
    os << "      </xs:extension>" << endl;
    os << "    </xs:complexContent>" << endl;
    os << "  </xs:complexType>" << endl;
    os << "  <xs:element name=\"multi_polygon\" "
                 "type=\"ogr:multi_polygon_Type\" "
                 "substitutionGroup=\"gml:_Feature\"/>" << endl;
    os << "  <xs:complexType name=\"multi_polygon_Type\">" << endl;
    os << "    <xs:complexContent>" << endl;
    os << "      <xs:extension base=\"gml:AbstractFeatureType\">" << endl;
    os << "        <xs:sequence>" << endl;
    os << "          <xs:element name=\"geometryProperty\" "
                         "type=\"gml:MultiPolygonPropertyType\" "
                         "nillable=\"true\" minOccurs=\"0\" maxOccurs=\"1\"/>" << endl;
    os << "          <xs:element name=\"fid\" nillable=\"true\" minOccurs=\"0\" "
                         "maxOccurs=\"1\">" << endl;
    os << "            <xs:simpleType>" << endl;
    os << "              <xs:restriction base=\"xs:string\">" << endl;
    os << "              </xs:restriction>" << endl;
    os << "            </xs:simpleType>" << endl;
    os << "          </xs:element>" << endl;
    os << "          <xs:element name=\"ID\" nillable=\"true\" minOccurs=\"0\" "
                         "maxOccurs=\"1\">" << endl;
    os << "            <xs:simpleType>" << endl;
    os << "              <xs:restriction base=\"xs:integer\">" << endl;
    os << "                <xs:totalDigits value=\"16\"/>" << endl;
    os << "              </xs:restriction>" << endl;
    os << "            </xs:simpleType>" << endl;
    os << "          </xs:element>" << endl;
    os << "        </xs:sequence>" << endl;
    os << "      </xs:extension>" << endl;
    os << "    </xs:complexContent>" << endl;
    os << "  </xs:complexType>" << endl;
    os << "</xs:schema>";
 
    return os.str().c_str();
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::To180(geos::geom::MultiPolygon *poly360) {
    try {
      // Let's take the 360 pieces that sit between 180 and 360 and move them
      //   to -180 to 0. To accomplish this, make a poly that fits 180 -> 360
      //   degrees longitude and intersect (this gives us the pieces that sit
      //   >180). Move this intersection to the left. Then make a poly that fits
      //   0 to 180 and intersect with the original. These two combined are the
      //   result.
      geos::geom::CoordinateSequence leftOf180Pts;
      leftOf180Pts.add(geos::geom::Coordinate(0, -90));
      leftOf180Pts.add(geos::geom::Coordinate(0, 90));
      leftOf180Pts.add(geos::geom::Coordinate(180, 90));
      leftOf180Pts.add(geos::geom::Coordinate(180, -90));
      leftOf180Pts.add(geos::geom::Coordinate(0, -90));
 
      unique_ptr<geos::geom::LinearRing> leftOf180Geom =
          globalFactory->createLinearRing(leftOf180Pts);
 
      geos::geom::Polygon *leftOf180Poly =
          globalFactory->createPolygon(std::move(leftOf180Geom)).release();
 
      geos::geom::CoordinateSequence rightOf180Pts;
      rightOf180Pts.add(geos::geom::Coordinate(180, -90));
      rightOf180Pts.add(geos::geom::Coordinate(180, 90));
      rightOf180Pts.add(geos::geom::Coordinate(360, 90));
      rightOf180Pts.add(geos::geom::Coordinate(360, -90));
      rightOf180Pts.add(geos::geom::Coordinate(180, -90));
 
      unique_ptr<geos::geom::LinearRing> rightOf180Geom =
          globalFactory->createLinearRing(rightOf180Pts);
 
      geos::geom::Polygon *rightOf180Poly =
          globalFactory->createPolygon(std::move(rightOf180Geom)).release();
 
      geos::geom::Geometry *preserved = Intersect(leftOf180Poly, poly360);
      geos::geom::Geometry *moving = Intersect(rightOf180Poly, poly360);
 
      geos::geom::CoordinateSequence *movingPts = moving->getCoordinates().release();
      geos::geom::CoordinateSequence movedPts;
 
      for(unsigned int i = 0; i < movingPts->getSize(); i ++) {
        movedPts.add(geos::geom::Coordinate(movingPts->getAt(i).x - 360.0,
                                            movingPts->getAt(i).y));
      }
 
      if(movedPts.getSize()) {
        movedPts.add(geos::geom::Coordinate(movedPts.getAt(0).x,
                                             movedPts.getAt(0).y));
      }
 
      geos::geom::Geometry *moved = globalFactory->createPolygon(
          globalFactory->createLinearRing(movedPts)).release();
 
      std::vector<const geos::geom::Geometry *> geomsForCollection;
      geomsForCollection.push_back(preserved);
      geomsForCollection.push_back(moved);
 
      geos::geom::GeometryCollection *the180Polys =
          Isis::globalFactory->createGeometryCollection(geomsForCollection).release();
 
      geos::geom::MultiPolygon *result = MakeMultiPolygon(the180Polys);
      delete the180Polys;
      the180Polys = NULL;
 
      geos::geom::MultiPolygon *fixedResult = FixSeam(result);
      delete result;
      result = NULL;
 
      return fixedResult;
    }
    catch(geos::util::GEOSException *exc) {
      IString msg = "Conversion to 180 failed. The reason given was [" +
          IString(exc->what()) + "]";
      delete exc;
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    catch(...) {
      IString msg = "Conversion to 180 failed. Could not determine the reason";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }
 
 
  double PolygonTools::Thickness(const geos::geom::MultiPolygon *mpolygon) {
    const geos::geom::Envelope *envelope = mpolygon->getEnvelopeInternal();
 
    double x = fabs(envelope->getMaxX() - envelope->getMinX());
    double y = fabs(envelope->getMaxY() - envelope->getMinY());
    double extent = max(x, y);
 
    return mpolygon->getArea() / (extent * extent);
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::Geometry *geom) {
    geos::geom::MultiPolygon *spiked = MakeMultiPolygon(geom);
    geos::geom::MultiPolygon *despiked = Despike(spiked);
 
    delete spiked;
    spiked = NULL;
 
    return despiked;
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::Despike(const geos::geom::MultiPolygon *multiPoly) {
    // Despike each polygon in the multipolygon
    vector<const geos::geom::Geometry *> *newPolys = new vector<const geos::geom::Geometry *>;
    for(unsigned int g = 0; g < multiPoly->getNumGeometries(); ++g) {
      const geos::geom::Polygon *poly =
          dynamic_cast<const geos::geom::Polygon *>(multiPoly->getGeometryN(g));
 
      // Despike each hole inside this polygon
      vector<unique_ptr<geos::geom::LinearRing>> holes;
      for(unsigned int h = 0; h < poly->getNumInteriorRing(); ++h) {
        const geos::geom::LineString *ls = poly->getInteriorRingN(h);
        unique_ptr<geos::geom::LinearRing> lr;
 
        // If the hole is not valid fix it
        // If the hole is NOT valid despike it
        lr = Despike(ls)->clone();
 
        if(!lr->isValid()) {
          geos::geom::LinearRing *fixed = FixGeometry(lr.release());
          lr.reset(fixed->clone().release());
        }
 
        // Save this hole if it is not too small
        if(!lr->isEmpty()) {
          holes.push_back(lr->clone());
        }
        else {
          lr.reset();
        }
      } // End holes loop
 
      // Work on the main polygon
      unique_ptr<geos::geom::LinearRing> ls = poly->getExteriorRing()->clone();
      unique_ptr<geos::geom::LinearRing> lr;
 
      lr = Despike(ls.release())->clone();
 
      try {
        if(!lr->isValid()) {
          geos::geom::LinearRing *fixed = FixGeometry(lr.release());
          lr.reset(fixed->clone().release());
        }
      }
      catch(IException &e) {
        // Sometimes despike and fix fail, but the input is really valid. We can just go
        // with the non-despiked polygon.
        if(ls->isValid() && ls->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
          lr.reset(ls.release());
        }
        else {
          throw;
        }
      }
 
      // Create a new polygon with the holes and save it
      if(!lr->isEmpty()) {
        geos::geom::Polygon *tp = Isis::globalFactory->createPolygon(std::move(lr), std::move(holes)).release();
 
        if(tp->isEmpty() || !tp->isValid()) {
          delete tp;
          newPolys->push_back(poly->clone().release());
        }
        else {
          newPolys->push_back(tp);
        }
      }
    } // End polygons loop
 
    // Create a new multipoly from the polygon(s)
    geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(*newPolys).release();
 
    if(!mp->isValid() || mp->isEmpty()) {
      delete mp;
      mp = NULL;
      IString msg = "Despike failed to correct the polygon";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    // if multipoly area changes more than 25% we did something bad to the multipolygon
    if(fabs((mp->getArea() / multiPoly->getArea()) - 1.0) > 0.50) {
      IString msg = "Despike failed to correct the polygon " + mp->toString();
      delete mp;
 
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    return mp;
  }
 
 
  geos::geom::LinearRing *PolygonTools::Despike(const geos::geom::LineString *lineString) {
    geos::geom::CoordinateSequence *vertices =
        new geos::geom::CoordinateSequence(*(lineString->getCoordinates()));
 
    // We need a full polygon to despike = 3 points (+1 for end==first) = at least 4 points
    if(vertices->getSize() < 4) {
      delete vertices;
      return Isis::globalFactory->createLinearRing(geos::geom::CoordinateSequence()).release();
    }
 
    // delete one of the duplicate first/end coordinates,
    //   spikes can occur here and the duplicate points throws off the test
    // deleteAt was removed in geos 3.8, this fix needs a better answer
    {
      std::vector<geos::geom::Coordinate> coords;
      vertices->toVector(coords);
      coords.pop_back();
      vertices->setPoints(coords);
      //vertices->deleteAt(vertices->getSize() - 1);
    }
 
    // Index will become negative if our first points are spiked, so we need
    //   to make it big enough to encapsulate the (unsigned int) from getSize()
    //   and allow negative numbers.
    for(long index = 0; index < (long)vertices->getSize(); index++) {
      // If we're under 3 vertices, we've despiked all the points out of the polygon :(
      if(vertices->getSize() < 3) {
        break;
      }
 
      // These are the array indices that we're going to test for a spike,
      //   the middle one is the one that's spiked if IsSpiked(...) is true.
      long testCoords[3] = {
        index - 1,
        index,
        index + 1
      };
 
      // Make sure the index is inside of our coordinates array (we'll have both too small/too
      // big of an index)
      for(int j = 0; j < 3; j++) {
        while(testCoords[j] < 0) {
          testCoords[j] += vertices->getSize();
        }
 
        while(testCoords[j] >= (long)vertices->getSize()) {
          testCoords[j] -= vertices->getSize();
        }
      }
 
      // Test the middle point for a spike
      if(IsSpiked(vertices->getAt(testCoords[0]),
                  vertices->getAt(testCoords[1]),
                  vertices->getAt(testCoords[2]))) {
        // It's spiked, delete it
        std::vector<geos::geom::Coordinate> coords;
        vertices->toVector(coords);
        coords.erase(coords.begin()+testCoords[1]);
        vertices->setPoints(coords);
 
        // Back up to the first test that is affected by this change
        index -= 2;
      }
    }
 
    if(vertices->getSize() < 3) {
      delete vertices;
      vertices = NULL;
 
      return Isis::globalFactory->createLinearRing(geos::geom::CoordinateSequence()).release();
    }
    else {
      // Duplicate the first vertex as the last to close the polygon
      vertices->add(vertices->getAt(0));
      return Isis::globalFactory->createLinearRing(*vertices).release();
    }
  }
 
  bool PolygonTools::IsSpiked(geos::geom::Coordinate first,
                              geos::geom::Coordinate middle,
                              geos::geom::Coordinate last) {
    return TestSpiked(first, middle, last) || TestSpiked(last, middle, first);
  }
 
  bool PolygonTools::TestSpiked(geos::geom::Coordinate first, geos::geom::Coordinate middle,
                                geos::geom::Coordinate last) {
    geos::geom::Point *middlePt = Isis::globalFactory->createPoint(middle).release();
    geos::geom::Point *lastPt = Isis::globalFactory->createPoint(last).release();
 
    geos::geom::CoordinateSequence coords;
    coords.add(first);
    coords.add(middle);
    geos::geom::LineString *line = Isis::globalFactory->createLineString(coords).release(); // line takes ownership
 
    // The lower the tolerance, the less this algorithm removes and thus
    //   the better chance of success in findimageoverlaps. However, if you
    //   lower the tolerance then there is also a greater chance of programs
    //   such as autoseed failing. 1% is the current tolerance.
    double tolerance = line->getLength() / 100.0;
 
    bool spiked = true;
 
    double distanceLastMiddle = geos::operation::distance::DistanceOp::distance(lastPt, middlePt);
    double distanceLastLine = geos::operation::distance::DistanceOp::distance(lastPt, line);
 
    if(distanceLastMiddle == 0.0) return true; // get rid of same point
 
    // Checks the ratio of the distance between the last point and the line, and the last point
    // and the middle point if the ratio is very small then there is a spike
    if(distanceLastLine / distanceLastMiddle >= .05) {
      spiked = false;
    }
 
    // If the point is away from the line, keep it
    if(spiked && distanceLastLine > tolerance) {
      spiked = false;
    }
 
    if(!spiked) {
      geos::geom::CoordinateSequence coords;
      coords.add(first);
      coords.add(middle);
      coords.add(last);
      coords.add(first);
 
      // shell takes ownership of coords
      geos::geom::LinearRing *shell = Isis::globalFactory->createLinearRing(coords).release();
      std::vector<geos::geom::LinearRing *> *empty = new std::vector<geos::geom::LinearRing *>;
 
      // poly takes ownership of shell and empty
      geos::geom::Polygon *poly = Isis::globalFactory->createPolygon(*shell, *empty);
 
 
      // if these 3 points define a straight line then the middle is worthless (defines nothing)
      // or problematic
      if(poly->getArea() < 1.0e-10) {
        spiked = true;
      }
 
      delete poly;
    }
 
 
    delete middlePt;
    delete lastPt;
    delete line;
 
    return spiked;
  }
 
 
  geos::geom::Geometry *PolygonTools::Intersect(const geos::geom::Geometry *geom1,
                                                const geos::geom::Geometry *geom2) {
    try {
      return Operate(geom1, geom2, (unsigned int)geos::operation::overlayng::OverlayNG::INTERSECTION);
    }
    catch(geos::util::GEOSException *exc) {
      IString msg = "Intersect operation failed. The reason given was [" + IString(exc->what()) + "]";
      delete exc;
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    catch(IException &e) {
      IString msg = "Intersect operation failed";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
    catch(...) {
      IString msg = "Intersect operation failed for an unknown reason";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }
 
 
  geos::geom::Geometry *PolygonTools::Operate(const geos::geom::Geometry *geom1,
                                              const geos::geom::Geometry *geom2,
                                              unsigned int opcode) {
 
    // geos::operation::overlayng::OverlayNG::OpCode code =
    //   (geos::operation::overlayng::OverlayNG::OpCode)opcode;
 
    geos::geom::Geometry *result = NULL;
    bool failed = true;
    geos::geom::Geometry *geomFirst  = MakeMultiPolygon(geom1);
    geos::geom::Geometry *geomSecond = MakeMultiPolygon(geom2);
 
    geos::operation::overlay::snap::GeometrySnapper snap(*geomFirst);
    geos::geom::Geometry *geomSnapped = snap.snapTo(*geomSecond, 1.0e-10)->clone().release();
    if(!geomSnapped->isValid()) {
      delete geomSnapped;
    }
    else {
      delete geomFirst;
      geomFirst = geomSnapped;
    }
 
    unsigned int precision = 15;
    unsigned int minPrecision = 13;
    while(failed) {
      try {
        // C++11: the geos BinaryOp returns an auto_ptr, we use release() to create a unique_ptr.
        std::unique_ptr< geos::geom::Geometry > resultAuto(
          geos::operation::overlay::snap::SnapOverlayOp::overlayOp(*geomFirst, *geomSecond, opcode));
        failed = false;
        result = resultAuto->clone().release();
      }
      catch(geos::util::GEOSException *exc) {
        // Just in case the clone failed....
        if(!failed || precision == minPrecision) throw;
 
        delete exc;
      }
      catch(...) {
        if(precision == minPrecision) {
          IString msg = "An unknown geos error occurred";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
 
        if(!failed) {
          IString msg = "An unknown geos error occurred when attempting to clone a geometry";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
      }
 
      // try reducing the precision
      if(failed) {
        precision --;
        geos::geom::Geometry *tmp = ReducePrecision(geomFirst, precision);
        delete geomFirst;
        geomFirst = tmp;
 
        tmp = ReducePrecision(geomSecond, precision);
        delete geomSecond;
        geomSecond = tmp;
      }
    }
 
    delete geomFirst;
    geomFirst = NULL;
 
    delete geomSecond;
    geomSecond = NULL;
 
    if(result && !result->isValid()) {
      try {
        geos::geom::Geometry *newResult = FixGeometry(result);
 
        if(fabs(newResult->getArea() / result->getArea() - 1.0) > 0.50) {
          delete newResult;
          delete result;
 
          IString msg = "Operation [" + IString((int)opcode) + "] failed";
          throw IException(IException::Programmer, msg, _FILEINFO_);
        }
 
        delete result;
        result = newResult;
      }
      catch(IException &e) {
        IString msg = "Operation [" + IString((int)opcode) + "] failed";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
    }
 
    if(result == NULL) {
      IString msg = "Operation [" + IString((int)opcode) + " failed";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    return result;
  }
 
  geos::geom::Geometry *PolygonTools::FixGeometry(const geos::geom::Geometry *geom) {
    if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
      return FixGeometry(dynamic_cast<const geos::geom::MultiPolygon *>(geom));
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
      return FixGeometry(dynamic_cast<const geos::geom::LinearRing *>(geom));
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
      return FixGeometry(dynamic_cast<const geos::geom::Polygon *>(geom));
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
      return FixGeometry(MakeMultiPolygon(geom));
    }
    else {
      IString msg = "PolygonTools::FixGeometry does not support [" + GetGeometryName(geom) + "]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }
 
  geos::geom::MultiPolygon *PolygonTools::FixGeometry(const geos::geom::MultiPolygon *poly) {
    // Maybe two points are on top of each other
    vector<const geos::geom::Geometry *> *newPolys = new vector<const geos::geom::Geometry *>;
 
    // Convert each polygon in this multi-polygon
    for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
      geos::geom::Polygon *fixedpoly = FixGeometry(
          dynamic_cast<const geos::geom::Polygon *>(
            poly->getGeometryN(geomIndex)));
      if(fixedpoly->isValid()) {
        newPolys->push_back(fixedpoly);
      }
      else {
        delete fixedpoly;
      }
      fixedpoly = NULL;
    }
 
    geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(*newPolys).release();
    return mp;
  }
 
 
  geos::geom::Polygon *PolygonTools::FixGeometry(const geos::geom::Polygon *poly) {
 
    // Convert each hole inside this polygon
    vector<unique_ptr<geos::geom::LinearRing>> holes;
    for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
      const geos::geom::LinearRing *thisHole = poly->getInteriorRingN(holeIndex);
 
      // We hope they are all linear rings (closed/simple), but if not just leave it be
      if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
        holes.push_back(thisHole->clone());
//      newRing = dynamic_cast<geos::geom::LinearRing *>(ring->clone());
        continue;
      }
 
      try {
        geos::geom::LinearRing *newHole = FixGeometry((geos::geom::LinearRing *)thisHole);
        holes.push_back(newHole->clone());
      }
      catch (IException &e) {
        IString msg = "Failed when attempting to fix interior ring of multipolygon";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
    } // end num holes in polygon loop
 
 
    const geos::geom::LineString *exterior = poly->getExteriorRing();
 
    try {
      unique_ptr<geos::geom::LinearRing> newExterior;
 
      if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
        newExterior.reset(FixGeometry((geos::geom::LinearRing *)exterior));
      }
      else {
        IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
                      "ring is not simple and closed";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
 
      return globalFactory->createPolygon(std::move(newExterior), std::move(holes)).release();
    }
    catch (IException &e) {
      IString msg = "Failed when attempting to fix exterior ring of polygon";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
  }
 
 
  geos::geom::LinearRing *PolygonTools::FixGeometry(const geos::geom::LinearRing *ring) {
 
    geos::geom::CoordinateSequence *coords = ring->getCoordinates().release();
 
    // Check this, just in case
    if(coords->getSize() < 4) {
      return globalFactory->createLinearRing(geos::geom::CoordinateSequence()).release();
    }
 
    geos::geom::CoordinateSequence *newCoords = new geos::geom::CoordinateSequence();
    const geos::geom::Coordinate *lastCoordinate = &coords->getAt(0);
    newCoords->add(*lastCoordinate);
 
    // Convert each coordinate in this hole
    for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
      const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
 
      // we're going to compare the decimal place of the current point to the decimal place
      //   of the difference, if they are drastically different then geos might not be seeing them
      //   correctly.
      double difference[2] = {
        lastCoordinate->x - thisCoordinate->x,
        lastCoordinate->y - thisCoordinate->y,
      };
 
      // geos isnt differentiating between points this close
      double minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));
 
      minDiff = min(minDiff, fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1])));
 
      // Cases where the difference in one direction is exactly zero, and the other direction is
      // next to zero appear often enough (especially in despike).
      if(difference[0] == 0.0 && difference[1] != 0.0) {
        // subtracted the two points, got deltaX = 0.0, use the y difference decimal place
        minDiff = fabs(DecimalPlace(thisCoordinate->y) - DecimalPlace(difference[1]));
      }
      else if(difference[1] == 0.0 && difference[0] != 0.0) {
        // subtracted the two points, got deltaY = 0.0, use the x difference decimal place
        minDiff = fabs(DecimalPlace(thisCoordinate->x) - DecimalPlace(difference[0]));
      }
      else if(difference[0] == 0.0 && difference[1] == 0.0) {
        // subtracted the two points, got 0.0, so it's same point... make sure it gets ignored!
        minDiff = 1E99;
      }
 
      // geos has a hard time differentiating when points get too close...
      if(minDiff <= 15) {
        newCoords->add(*thisCoordinate);
        lastCoordinate = thisCoordinate;
      }
    } // end num coords in hole loop
 
    newCoords->add(geos::geom::Coordinate(newCoords->getAt(0).x, newCoords->getAt(0).y));
    geos::geom::LinearRing *newRing = NULL;
 
    // Now that we've weeded out any bad coordinates, let's rebuild the geometry
    try {
      if(newCoords->getSize() > 3) {
        newRing = globalFactory->createLinearRing(*newCoords).release();
      }
      else {
        delete newCoords;
        newCoords = NULL;
      }
    }
    catch(geos::util::GEOSException *exc) {
      delete exc;
      exc = NULL;
 
      IString msg = "Error when attempting to fix linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if(newRing && !newRing->isValid() && ring->isValid()) {
      IString msg = "Failed when attempting to fix linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    else if(!newRing || !newRing->isValid()) {
      if(newRing) {
        delete newRing;
      }
 
      newRing = dynamic_cast<geos::geom::LinearRing *>(ring->clone().release());
    }
 
    return newRing;
  }
 
 
  int PolygonTools::DecimalPlace(double num) {
    // 0.0 = decimal place 0
    if(num == 0.0) return 0;
 
    num = fabs(num);
 
    int decimalPlace = 1;
    while(num < 1.0) {
      num *= 10.0;
      decimalPlace --;
    }
 
    while(num > 10.0) {
      num /= 10.0;
      decimalPlace ++;
    }
 
    return decimalPlace;
  }
 
  geos::geom::Geometry *PolygonTools::Difference(const geos::geom::Geometry *geom1,
                                                 const geos::geom::Geometry *geom2) {
    try {
      return Operate(geom1, geom2, (unsigned int)geos::operation::overlayng::OverlayNG::DIFFERENCE);
    }
    catch(geos::util::GEOSException *exc) {
      IString msg = "Difference operation failed. The reason given was [" +
                    IString(exc->what()) + "]";
      delete exc;
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    catch(IException &e) {
      IString msg = "Difference operation failed";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
    catch(...) {
      IString msg = "Difference operation failed for an unknown reason";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::MakeMultiPolygon(const geos::geom::Geometry *geom) {
    // The area of the geometry is too small, so just ignore it.
    if(geom->isEmpty()) {
      return Isis::globalFactory->createMultiPolygon().release();
    }
 
    else if(geom->getArea() - DBL_EPSILON <= DBL_EPSILON) {
      return Isis::globalFactory->createMultiPolygon().release();
    }
 
    else if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
      return dynamic_cast<geos::geom::MultiPolygon *>(geom->clone().release());
    }
 
    else if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
      vector<const geos::geom::Geometry *> polys;
      polys.push_back(geom->clone().release());
      geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(polys).release();
      return mp;
    }
 
    else if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
      vector<const geos::geom::Geometry *> polys;
      const geos::geom::GeometryCollection *gc =
          dynamic_cast<const geos::geom::GeometryCollection *>(geom);
      for(unsigned int i = 0; i < gc->getNumGeometries(); ++i) {
        geos::geom::MultiPolygon *subMultiPoly =
            MakeMultiPolygon(gc->getGeometryN(i));
 
        for(unsigned int subPoly = 0;
            subPoly < subMultiPoly->getNumGeometries();
            subPoly ++) {
          const geos::geom::Polygon *poly =
              dynamic_cast<const geos::geom::Polygon *>(
                subMultiPoly->getGeometryN(subPoly));
          polys.push_back(poly->clone().release());
        }
      }
 
      geos::geom::MultiPolygon *mp =
          Isis::globalFactory->createMultiPolygon(polys).release();
      if(mp->getArea() - DBL_EPSILON <= DBL_EPSILON) {
        delete mp;
        mp = Isis::globalFactory->createMultiPolygon().release();
      }
 
      return mp;
    }
 
    // All other geometry types are invalid so ignore them
    else {
      return Isis::globalFactory->createMultiPolygon().release();
    }
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::FixSeam(
      const geos::geom::Polygon *polyA, const geos::geom::Polygon *polyB) {
    geos::geom::CoordinateSequence *polyAPoints = polyA->getCoordinates().release();
    geos::geom::CoordinateSequence *polyBPoints = polyB->getCoordinates().release();
 
    unsigned int aIntersectionBegin = 0;
    unsigned int aIntersectionEnd = 0;
    unsigned int bIntersectionBegin = 0;
    unsigned int bIntersectionEnd = 0;
 
    bool intersectionStarted = false;
    bool intersectionEnded = false;
 
    unsigned int lastBMatch  = 0;
    for (unsigned int i = 0;
        !intersectionEnded && i < polyAPoints->getSize();
        i++) {
 
      bool foundEquivalent = false;
 
      geos::geom::Coordinate coordA = polyAPoints->getAt(i);
      coordA = *ReducePrecision(&coordA, 13);
 
      for (unsigned int j = 0;
           !foundEquivalent && j < polyBPoints->getSize();
           j++) {
        geos::geom::Coordinate coordB = polyBPoints->getAt(j);
        coordB = *ReducePrecision(&coordB, 13);
 
        foundEquivalent = coordA.equals(coordB);
 
        if (foundEquivalent) lastBMatch = j;
 
        if (foundEquivalent && !intersectionStarted) {
          intersectionStarted = true;
          aIntersectionBegin = i;
          bIntersectionBegin = j;
        }
      }
 
      if (!foundEquivalent && intersectionStarted && !intersectionEnded) {
        intersectionEnded = true;
        aIntersectionEnd = i;
        bIntersectionEnd = lastBMatch;
      }
    }
 
    geos::geom::MultiPolygon * result = NULL;
    if (intersectionStarted && intersectionEnded) {
      geos::geom::CoordinateSequence merged;
 
      unsigned int i = 0;
      for (i = 0; i < aIntersectionBegin; i ++) {
        merged.add(polyAPoints->getAt(i));
      }
 
      i = bIntersectionBegin;
      while (i != bIntersectionEnd) {
        merged.add(polyBPoints->getAt(i));
        i++;
        if (i >= polyBPoints->getSize()) i = 0;
      }
 
      for (i = aIntersectionEnd; i < polyAPoints->getSize() - 1; i++) {
        merged.add(polyAPoints->getAt(i));
      }
 
      merged.add(merged.getAt(0));
      result = MakeMultiPolygon(globalFactory->createPolygon(
          globalFactory->createLinearRing(merged)).release());
    }
 
    return result;
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::FixSeam(
      const geos::geom::MultiPolygon *poly) {
 
    std::vector<const geos::geom::Geometry *> polys;
 
 
    for(unsigned int copyIndex = 0;
        copyIndex < poly->getNumGeometries();
        copyIndex ++) {
      polys.push_back(poly->getGeometryN(copyIndex)->clone().release());
    }
 
    unsigned int outerPolyIndex = 0;
 
    while(outerPolyIndex + 1 < polys.size()) {
      unsigned int innerPolyIndex = outerPolyIndex + 1;
 
      while(innerPolyIndex < polys.size()) {
        geos::geom::MultiPolygon *fixedPair = FixSeam(
            dynamic_cast<const geos::geom::Polygon *>(polys.at(outerPolyIndex)),
            dynamic_cast<const geos::geom::Polygon *>(polys.at(innerPolyIndex)));
 
        if(fixedPair != NULL) {
          const geos::geom::Geometry *oldInnerPoly = polys.at(innerPolyIndex);
          const geos::geom::Geometry *oldOuterPoly = polys.at(outerPolyIndex);
 
          polys.erase(polys.begin() + innerPolyIndex);
          polys[outerPolyIndex] = fixedPair->getGeometryN(0)->clone().release();
          innerPolyIndex = outerPolyIndex + 1;
 
          delete fixedPair;
          fixedPair = NULL;
 
          delete oldInnerPoly;
          oldInnerPoly = NULL;
 
          delete oldOuterPoly;
          oldOuterPoly = NULL;
        }
        else {
          innerPolyIndex ++;
        }
      }
 
      outerPolyIndex ++;
    }
 
    return globalFactory->createMultiPolygon(polys).release();
  }
 
 
  geos::geom::Geometry *PolygonTools::ReducePrecision(const geos::geom::Geometry *geom,
                                                      unsigned int precision) {
    if(geom->getGeometryTypeId() == geos::geom::GEOS_MULTIPOLYGON) {
      return ReducePrecision(
          dynamic_cast<const geos::geom::MultiPolygon *>(geom), precision);
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
      return ReducePrecision(
          dynamic_cast<const geos::geom::LinearRing *>(geom), precision);
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_POLYGON) {
      return ReducePrecision(
          dynamic_cast<const geos::geom::Polygon *>(geom), precision);
    }
    if(geom->getGeometryTypeId() == geos::geom::GEOS_GEOMETRYCOLLECTION) {
      return ReducePrecision(MakeMultiPolygon(geom), precision);
    }
    else {
      IString msg = "PolygonTools::ReducePrecision does not support [" +
                    GetGeometryName(geom) + "]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::ReducePrecision(const geos::geom::MultiPolygon *poly,
                                                          unsigned int precision) {
    // Maybe two points are on top of each other
    vector<const geos::geom::Geometry *> newPolys;
 
    // Convert each polygon in this multi-polygon
    for(unsigned int geomIndex = 0; geomIndex < poly->getNumGeometries(); geomIndex ++) {
      geos::geom::Geometry *lowerPrecision = ReducePrecision(
          dynamic_cast<const geos::geom::Polygon *>(
            poly->getGeometryN(geomIndex)),
          precision);
 
      if(!lowerPrecision->isEmpty()) {
        newPolys.push_back(lowerPrecision);
      }
      else {
        delete lowerPrecision;
      }
    }
 
    geos::geom::MultiPolygon *mp = Isis::globalFactory->createMultiPolygon(newPolys).release();
    return mp;
  }
 
 
  geos::geom::Polygon *PolygonTools::ReducePrecision(const geos::geom::Polygon *poly,
                                                     unsigned int precision) {
    // Convert each hole inside this polygon
    vector<unique_ptr<geos::geom::LinearRing>> holes;
    for(unsigned int holeIndex = 0; holeIndex < poly->getNumInteriorRing(); holeIndex ++) {
      const geos::geom::LinearRing *thisHole = poly->getInteriorRingN(holeIndex);
 
      // We hope they are all linear rings (closed/simple), but if not just leave it be
      if(thisHole->getGeometryTypeId() != geos::geom::GEOS_LINEARRING) {
        holes.push_back(thisHole->clone());
        continue;
      }
 
      try {
        geos::geom::LinearRing *newHole = ReducePrecision((geos::geom::LinearRing *)thisHole,
                                                          precision);
 
        if(!newHole->isEmpty()) {
          holes.push_back(newHole->clone());
        }
        else {
          delete newHole;
        }
 
      }
      catch(IException &e) {
        IString msg = "Failed when attempting to fix interior ring of multipolygon";
        throw IException(e, IException::Programmer, msg, _FILEINFO_);
      }
    } // end num holes in polygon loop
 
 
    const geos::geom::LineString *exterior = poly->getExteriorRing();
 
    try {
      unique_ptr<geos::geom::LinearRing> newExterior;
 
      if(exterior->getGeometryTypeId() == geos::geom::GEOS_LINEARRING) {
        newExterior.reset(ReducePrecision((geos::geom::LinearRing *)exterior, precision));
      }
      else {
        IString msg = "Failed when attempting to fix exterior ring of polygon. The exterior "
                      "ring is not simple and closed";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
 
      return globalFactory->createPolygon(std::move(newExterior), std::move(holes)).release();
    }
    catch(IException &e) {
      IString msg = "Failed when attempting to fix exterior ring of polygon";
      throw IException(e, IException::Programmer, msg, _FILEINFO_);
    }
  }
 
 
  geos::geom::LinearRing *PolygonTools::ReducePrecision(const geos::geom::LinearRing *ring,
                                                        unsigned int precision) {
    geos::geom::CoordinateSequence *coords = ring->getCoordinates().release();
 
    // Check this, just in case
    if(coords->getSize() <= 0) {
      return dynamic_cast<geos::geom::LinearRing *>(ring->clone().release());
    }
 
    geos::geom::CoordinateSequence newCoords;
    geos::geom::Coordinate *coord = ReducePrecision(&coords->getAt(0), precision);
    newCoords.add(*coord);
    delete coord;
    coord = NULL;
 
    // Convert each coordinate in this ring
    for(unsigned int coordIndex = 1; coordIndex < coords->getSize() - 1; coordIndex ++) {
      const geos::geom::Coordinate *thisCoordinate = &coords->getAt(coordIndex);
      coord = ReducePrecision(thisCoordinate, precision);
      newCoords.add(*coord);
      delete coord;
      coord = NULL;
    }
 
    newCoords.add(geos::geom::Coordinate(newCoords.getAt(0).x, newCoords.getAt(0).y));
    geos::geom::LinearRing *newRing = NULL;
 
    // Now that we've weeded out any bad coordinates, let's rebuild the geometry
    try {
      newRing = globalFactory->createLinearRing(newCoords).release();
    }
    catch(geos::util::GEOSException *exc) {
      delete exc;
      exc = NULL;
 
      IString msg = "Error when attempting to reduce precision of linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    // try to despike
    try {
      geos::geom::LinearRing *tmp = Despike(newRing);
      delete newRing;
      newRing = tmp;
    }
    catch(IException &e) {
    }
 
    if(!newRing->isValid()) {
      IString msg = "Failed when attempting to reduce precision of linear ring";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    return newRing;
  }
 
 
  geos::geom::Coordinate *PolygonTools::ReducePrecision(const geos::geom::Coordinate *coord,
                                                        unsigned int precision) {
    return new geos::geom::Coordinate(
             ReducePrecision(coord->x, precision),
             ReducePrecision(coord->y, precision),
             ReducePrecision(coord->z, precision)
           );
  }
 
 
  double PolygonTools::ReducePrecision(double num, unsigned int precision) {
    double result = num;
 
    // If num == nan then this test will fail
    if (num == num) {
      int decimalPlace = DecimalPlace(num);
      double factor = pow(10.0, (int)decimalPlace);
 
      // reduced num is in the form 0.nnnnnnnnnn...
      double reducedNum = num / factor;
 
      double cutoff = pow(10.0, (int)precision);
      double round_offset = (num < 0) ? -0.5 : 0.5;
 
      // cast off the digits past the precision's place
      reducedNum = ((long long)(reducedNum * cutoff + round_offset)) / cutoff;
      result = reducedNum * factor;
    }
 
    return result;
  }
 
 
  QString PolygonTools::GetGeometryName(const geos::geom::Geometry *geom) {
    switch(geom->getGeometryTypeId()) {
      case geos::geom::GEOS_POINT:
        return "Point";
      case geos::geom::GEOS_LINESTRING:
        return "Line String";
      case geos::geom::GEOS_LINEARRING:
        return "Linear Ring";
      case geos::geom::GEOS_POLYGON:
        return "Polygon";
      case geos::geom::GEOS_MULTIPOINT:
        return "Multi Point";
      case geos::geom::GEOS_MULTILINESTRING:
        return "Multi Line String";
      case geos::geom::GEOS_MULTIPOLYGON:
        return "Multi Polygon";
      case geos::geom::GEOS_GEOMETRYCOLLECTION:
        return "Geometry Collection";
      default:
        return "UNKNOWN";
    }
  }
 
 
  bool PolygonTools::Equal(const geos::geom::MultiPolygon *poly1,
                           const geos::geom::MultiPolygon *poly2) {
 
    vector<const geos::geom::Polygon *> polys1;
    vector<const geos::geom::Polygon *> polys2;
 
    if(poly1->getNumGeometries() != poly2->getNumGeometries())  return false;
 
    // Convert each polygon in this multi-polygon
    for(unsigned int geomIndex = 0; geomIndex < poly1->getNumGeometries(); geomIndex ++) {
      polys1.push_back(dynamic_cast<const geos::geom::Polygon *>(
          poly1->getGeometryN(geomIndex)));
      polys2.push_back(dynamic_cast<const geos::geom::Polygon *>(
          poly2->getGeometryN(geomIndex)));
    }
 
    for(int p1 = polys1.size() - 1; (p1 >= 0) && polys1.size(); p1 --) {
      for(int p2 = polys2.size() - 1; (p2 >= 0) && polys2.size(); p2 --) {
        if(Equal(polys1[p1], polys2[p2])) {
          // Delete polys1[p1] by replacing it with the last Polygon in polys1
          polys1[p1] = polys1[polys1.size()-1];
          polys1.resize(polys1.size() - 1);
          // Delete polys2[p2] by replacing it with the last Polygon in polys2
          polys2[p2] = polys2[polys2.size()-1];
          polys2.resize(polys2.size() - 1);
        }
      }
    }
 
    return (polys1.size() == 0) && (polys2.size() == 0);
  }
 
 
  bool PolygonTools::Equal(const geos::geom::Polygon *poly1, const geos::geom::Polygon *poly2) {
    vector<const geos::geom::LineString *> holes1;
    vector<const geos::geom::LineString *> holes2;
 
    if(poly1->getNumInteriorRing() != poly2->getNumInteriorRing())  return false;
 
    if(!Equal(poly1->getExteriorRing(), poly2->getExteriorRing()))  return false;
 
    // Convert each hole inside this polygon
    for(unsigned int holeIndex = 0; holeIndex < poly1->getNumInteriorRing(); holeIndex ++) {
 
      // We hope they are all linear rings (closed/simple), but if not just leave it be
      if(poly1->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
        holes1.push_back(poly1->getInteriorRingN(holeIndex));
      }
 
      if(poly2->getInteriorRingN(holeIndex)->getGeometryTypeId() == geos::geom::GEOS_LINESTRING) {
        holes2.push_back(poly2->getInteriorRingN(holeIndex));
      }
 
    }
 
    if(holes1.size() != holes2.size())  return false;
 
    for(int h1 = holes1.size() - 1; (h1 >= 0) && holes1.size(); h1 --) {
      for(int h2 = holes2.size() - 1; (h2 >= 0) && holes2.size(); h2 --) {
        if(Equal(holes1[h1], holes2[h2])) {
          // Delete holes1[h1] by replacing it with the last Polygon in holes1
          holes1[h1] = holes1[holes1.size()-1];
          holes1.resize(holes1.size() - 1);
          // Delete holes2[h2] by replacing it with the last Polygon in holes2
          holes2[h2] = holes2[holes2.size()-1];
          holes2.resize(holes2.size() - 1);
        }
      }
    }
 
    return (holes1.size() == 0) && (holes2.size() == 0);
  }
 
 
  bool PolygonTools::Equal(const geos::geom::LineString *lineString1,
                           const geos::geom::LineString *lineString2) {
 
    geos::geom::CoordinateSequence *coords1 = lineString1->getCoordinates().release();
    geos::geom::CoordinateSequence *coords2 = lineString2->getCoordinates().release();
    bool isEqual = true;
 
    if(coords1->getSize() != coords2->getSize()) isEqual = false;
 
    unsigned int index1 = 0;
    unsigned int index2 = 0;
 
    // -1 extra for dupicate start/end coordinates
    for(; index2 < coords2->getSize() - 1 && isEqual; index2 ++) {
      if(Equal(coords1->getAt(index1), coords2->getAt(index2)))  break;
    }
 
    if(index2 == coords2->getSize() - 1) isEqual = false;
 
    for(; index1 < coords1->getSize() - 1 && isEqual; index1 ++, index2 ++) {
      if(!Equal(coords1->getAt(index1), coords2->getAt(index2 % (coords2->getSize() - 1)))) {
        isEqual = false;
      }
    }
 
    delete coords1;
    delete coords2;
    return isEqual;
  }
 
 
  bool PolygonTools::Equal(const geos::geom::Coordinate &coord1,
                           const geos::geom::Coordinate &coord2) {
 
    if(!Equal(coord1.x, coord2.x))  return false;
    if(!Equal(coord1.y, coord2.y))  return false;
    if(!Equal(coord1.y, coord2.y))  return false;
 
    return true;
  }
 
 
  bool PolygonTools::Equal(const double d1, const double d2) {
    const double cutoff = 1e15;
 
    if(DecimalPlace(d1) != DecimalPlace(d2)) return false;
 
    int decimalPlace = DecimalPlace(d1);
    double factor = pow(10.0, (int)decimalPlace);
 
    // reduced num is in the form 0.nnnnnnnnnn...
    double reducedNum = d1 / factor;
 
    double round_offset = (d1 < 0) ? -0.5 : 0.5;
 
    // cast off the digits past the precision's place
    long long num1 = ((long long)(reducedNum * cutoff + round_offset));
 
    factor = pow(10.0, (int)decimalPlace);
 
    // reduced num is in the form 0.nnnnnnnnnn...
    reducedNum = d2 / factor;
 
    round_offset = (d2 < 0) ? -0.5 : 0.5;
 
    // cast off the digits past the precision's place
    long long num2 = ((long long)(reducedNum * cutoff + round_offset));
 
 
    return (num1 == num2);
  }
 
 
  geos::geom::MultiPolygon *PolygonTools::SplitPolygonOn360(const geos::geom::Polygon *inPoly) {
    bool convertLon = false;
    bool negAdjust = false;
    bool newCoords = false;  //  coordinates have been adjusted
    geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateSequence();
    double lon, lat;
    double lonOffset = 0;
    geos::geom::CoordinateSequence *inPolyCoords = inPoly->getCoordinates().release();
    double prevLon = inPolyCoords->getAt(0).x;
    double prevLat = inPolyCoords->getAt(0).y;
 
    newLonLatPts->add(geos::geom::Coordinate(prevLon, prevLat));
    double dist = 0.0;
    for (unsigned int i = 1; i < inPolyCoords->getSize(); i++) {
      lon = inPolyCoords->getAt(i).x;
      lat = inPolyCoords->getAt(i).y;
      // check to see if you just crossed the Meridian
      if (abs(lon - prevLon) > 180 && (prevLat != 90 && prevLat != -90)) {
        newCoords = true;
        // if you were already converting then stop (crossed Meridian even number of times)
        if (convertLon) {
          convertLon = false;
          lonOffset = 0;
        }
        else {   // Need to start converting again, deside how to adjust coordinates
          if ((lon - prevLon) > 0) {
            lonOffset = -360.;
            negAdjust = true;
          }
          else if ((lon - prevLon) < 0) {
            lonOffset = 360.;
            negAdjust = false;
          }
          convertLon = true;
        }
      }
 
      // Change to a minimum calculation
      if (newCoords  &&  dist == 0.0) {
        double longitude = (lon + lonOffset) - prevLon;
        double latitude = lat - prevLat;
        dist = std::sqrt((longitude * longitude) + (latitude * latitude));
      }
 
      // add coord
      newLonLatPts->add(geos::geom::Coordinate(lon + lonOffset, lat));
 
      // set current to old
      prevLon = lon;
      prevLat = lat;
    }
 
    delete inPolyCoords;
 
    // Nothing was done so return
    if (!newCoords) {
      geos::geom::Polygon *newPoly = globalFactory->createPolygon
                                     (globalFactory->createLinearRing(*newLonLatPts)).release();
      geos::geom::MultiPolygon *multi_polygon = PolygonTools::MakeMultiPolygon(newPoly);
      delete newLonLatPts;
      return multi_polygon;
    }
 
    // bisect into seperate polygons
    try {
      geos::geom::Polygon *newPoly = globalFactory->createPolygon
                                     (globalFactory->createLinearRing(*newLonLatPts)).release();
 
      geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateSequence();
      geos::geom::CoordinateSequence *pts2 = new geos::geom::CoordinateSequence();
 
      // Depending on direction of compensation bound accordingly
      //***************************************************
 
      // please verify correct if you change these values
      //***************************************************
      if (negAdjust) {
        pts->add(geos::geom::Coordinate(0., 90.));
        pts->add(geos::geom::Coordinate(-360., 90.));
        pts->add(geos::geom::Coordinate(-360., -90.));
        pts->add(geos::geom::Coordinate(0., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts->add(geos::geom::Coordinate(0.0, lat));
        }
        pts->add(geos::geom::Coordinate(0., 90.));
        pts2->add(geos::geom::Coordinate(0., 90.));
        pts2->add(geos::geom::Coordinate(360., 90.));
        pts2->add(geos::geom::Coordinate(360., -90.));
        pts2->add(geos::geom::Coordinate(0., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts2->add(geos::geom::Coordinate(0.0, lat));
        }
        pts2->add(geos::geom::Coordinate(0., 90.));
      }
      else {
        pts->add(geos::geom::Coordinate(360., 90.));
        pts->add(geos::geom::Coordinate(720., 90.));
        pts->add(geos::geom::Coordinate(720., -90.));
        pts->add(geos::geom::Coordinate(360., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts->add(geos::geom::Coordinate(360.0, lat));
        }
        pts->add(geos::geom::Coordinate(360., 90.));
        pts2->add(geos::geom::Coordinate(360., 90.));
        pts2->add(geos::geom::Coordinate(0., 90.));
        pts2->add(geos::geom::Coordinate(0., -90.));
        pts2->add(geos::geom::Coordinate(360., -90.));
        for (double lat = -90.0 + dist; lat < 90.0; lat += dist) {
          pts2->add(geos::geom::Coordinate(360.0, lat));
        }
        pts2->add(geos::geom::Coordinate(360., 90.));
      }
 
      geos::geom::Polygon *boundaryPoly = globalFactory->createPolygon
                                          (globalFactory->createLinearRing(*pts)).release();
      geos::geom::Polygon *boundaryPoly2 = globalFactory->createPolygon
                                           (globalFactory->createLinearRing(*pts2)).release();
      /*------------------------------------------------------------------------
      /  Intersecting the original polygon (converted coordinates) with the
      /  boundary polygons will create the multi polygons with the converted coordinates.
      /  These will need to be converted back to the original coordinates.
      /-----------------------------------------------------------------------*/
      geos::geom::Geometry *intersection = PolygonTools::Intersect(newPoly, boundaryPoly);
      geos::geom::MultiPolygon *convertPoly = PolygonTools::MakeMultiPolygon(intersection);
      delete intersection;
 
      intersection = PolygonTools::Intersect(newPoly, boundaryPoly2);
      geos::geom::MultiPolygon *convertPoly2 = PolygonTools::MakeMultiPolygon(intersection);
      delete intersection;
 
      /*------------------------------------------------------------------------
      / Adjust points created in the negative space or >360 space to be back in
      / the 0-360 world.  This will always only need to be done on convertPoly.
      / Then add geometries to finalpolys.
      /-----------------------------------------------------------------------*/
      vector<const geos::geom::Geometry *> finalpolys;
      geos::geom::Geometry *newGeom = NULL;
 
      for (unsigned int i = 0; i < convertPoly->getNumGeometries(); i++) {
        newGeom = (convertPoly->getGeometryN(i))->clone().release();
        geos::geom::CoordinateSequence *pts3 = convertPoly->getGeometryN(i)->getCoordinates().release();
        geos::geom::CoordinateSequence *newLonLatPts = new geos::geom::CoordinateSequence();
 
        // fix the points
        for (size_t k = 0; k < pts3->getSize() ; k++) {
          double lon = pts3->getAt(k).x;
          double lat = pts3->getAt(k).y;
          if (negAdjust) {
            lon = lon + 360;
          }
          else {
            lon = lon - 360;
          }
          newLonLatPts->add(geos::geom::Coordinate(lon, lat), k);
        }
        // Add the points to polys
        finalpolys.push_back(globalFactory->createPolygon
                              (globalFactory->createLinearRing(*newLonLatPts)).release());
      }
 
      // This loop is over polygons that will always be in 0-360 space no need to convert
      for (unsigned int i = 0; i < convertPoly2->getNumGeometries(); i++) {
        newGeom = (convertPoly2->getGeometryN(i))->clone().release();
        finalpolys.push_back(newGeom);
      }
 
      geos::geom::MultiPolygon *multi_polygon = globalFactory->createMultiPolygon(finalpolys).release();
 
      delete newLonLatPts;
      delete pts;
      delete pts2;
      return multi_polygon;
    }
    catch(geos::util::IllegalArgumentException *geosIll) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "geos illegal argument [" + IString(geosIll->what()) + "]";
      delete geosIll;
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
    catch(geos::util::GEOSException *geosExc) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "geos exception [" + IString(geosExc->what()) + "]";
      delete geosExc;
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
    catch(IException &e) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "isis operation exception [" + IString(e.what()) + "]";
      throw IException(e, IException::Unknown, msg, _FILEINFO_);
    }
    catch(...) {
      std::string msg = "Unable to split polygon on longitude boundry (SplitPolygonOn360) due to ";
      msg += "unknown exception";
      throw IException(IException::Unknown, msg, _FILEINFO_);
    }
  }
} // end namespace isis