UniversalGroundMap.cpp

#include "UniversalGroundMap.h"

#include <QPointF>

#include "Camera.h"
#include "CameraFactory.h"
#include "ImagePolygon.h"
#include "Latitude.h"
#include "Longitude.h"
#include "Latitude.h"
#include "Longitude.h"
#include "PolygonTools.h"
#include "Projection.h"
#include "RingPlaneProjection.h"
#include "TProjection.h"
#include "ProjectionFactory.h"
#include "SurfacePoint.h"
#include "Target.h"

namespace Isis {
  UniversalGroundMap::UniversalGroundMap(Cube &cube, CameraPriority priority) {
    p_camera = NULL;
    p_projection = NULL;

    Pvl &pvl = *cube.label();
    try {
      if(priority == CameraFirst)
        p_camera = CameraFactory::Create(cube);
      else
        p_projection = Isis::ProjectionFactory::CreateFromCube(pvl);
    }
    catch (IException &firstError) {
      p_camera = NULL;
      p_projection = NULL;

      try {
        if(priority == CameraFirst)
          p_projection = Isis::ProjectionFactory::CreateFromCube(pvl);
        else
          p_camera = CameraFactory::Create(cube);
      }
      catch (IException &secondError) {
        p_projection = NULL;
        QString msg = "Could not create camera or projection for [" +
                          cube.fileName() + "]";
        IException realError(IException::Unknown, msg, _FILEINFO_);
        realError.append(firstError);
        realError.append(secondError);
        throw realError;
      }
    }
  }

  void UniversalGroundMap::SetBand(const int band) {
    if (p_camera != NULL)
      p_camera->SetBand(band);
  }



  UniversalGroundMap::~UniversalGroundMap() {
    if (p_camera != NULL) {
      delete p_camera;
      p_camera = NULL;
    }

    if (p_projection != NULL) {
      delete p_projection;
      p_projection = NULL;
    }
  }

  bool UniversalGroundMap::SetUniversalGround(double lat, double lon) {
    if (p_camera != NULL) {
      if (p_camera->SetUniversalGround(lat, lon)) {  // This should work for rings (radius,azimuth)
        return p_camera->InCube();
      }
      else {
        return false;
      }
    }
    else {
      return p_projection->SetUniversalGround(lat, lon); // This should work for rings (radius,azimuth)
    }
  }


  bool UniversalGroundMap::SetGround(Latitude lat, Longitude lon) {
    if(p_camera != NULL) {
      if(p_camera->SetGround(lat, lon)) {  // This should work for rings (radius,azimuth)
        return p_camera->InCube();
      }
      else {
        return false;
      }
    }
    else {
      double universalLat = lat.degrees();
      double universalLon = lon.degrees();
      return p_projection->SetUniversalGround(universalLat, universalLon);  // This should work for rings (radius,azimuth)
    }
  }


  bool UniversalGroundMap::SetUnboundGround(Latitude lat, Longitude lon) {
    if(p_camera != NULL) {
      if(p_camera->SetGround(lat, lon)) {  // This should work for rings (radius,azimuth)
        return p_camera->InCube();
      }
      else {
        return false;
      }
    }
    else {
      double universalLat = lat.degrees();
      double universalLon = lon.degrees();
      return p_projection->SetUnboundUniversalGround(universalLat, universalLon);
    }
  }


  bool UniversalGroundMap::SetGround(const SurfacePoint &sp) {
    if (p_camera != NULL) {
      if (p_camera->SetGround(sp)) {
        return p_camera->InCube();
      }
      else {
        return false;
      }
    }
    else {
      return p_projection->SetUniversalGround(sp.GetLatitude().degrees(),
                                              sp.GetLongitude().degrees());   // This should work for rings (radius,azimuth)
    }
  }

  double UniversalGroundMap::Sample() const {
    if (p_camera != NULL) {
      return p_camera->Sample();
    }
    else {
      return p_projection->WorldX();
    }
  }

  double UniversalGroundMap::Line() const {
    if (p_camera != NULL) {
      return p_camera->Line();
    }
    else {
      return p_projection->WorldY();
    }
  }

  bool UniversalGroundMap::SetImage(double sample, double line) {
    if (p_camera != NULL) {
      return p_camera->SetImage(sample, line);
    }
    else {
      return p_projection->SetWorld(sample, line);
    }
  }

  double UniversalGroundMap::UniversalLatitude() const {
    if (p_camera != NULL) {
      return p_camera->UniversalLatitude();
    }
    else {
      //  Is this a triaxial projection or ring projection.  If ring return Radius as latitude
      Projection::ProjectionType projType = p_projection->projectionType();
      if (projType == Projection::Triaxial) {
        TProjection *tproj = (TProjection *) p_projection;
        return tproj->UniversalLatitude();
      }
      else {
        RingPlaneProjection *rproj = (RingPlaneProjection *) p_projection;
        return rproj->RingRadius();
      }
    }
  }

  double UniversalGroundMap::UniversalLongitude() const {
    if (p_camera != NULL) {
      return p_camera->UniversalLongitude();
    }
    else {
      //  Is this a triaxial projection or ring projection.  If ring return ring longitude as
      //    longitude
      Projection::ProjectionType projType = p_projection->projectionType();
      if (projType == Projection::Triaxial) {
        TProjection *tproj = (TProjection *) p_projection;
        return tproj->UniversalLongitude();
      }
      else {
        RingPlaneProjection *rproj = (RingPlaneProjection *) p_projection;
        return rproj->RingLongitude();
      }
    }
  }

  double UniversalGroundMap::Resolution() const {
    if (p_camera != NULL) {
      return p_camera->PixelResolution();
    }
    else {
      return p_projection->Resolution();
    }
  }


  bool UniversalGroundMap::GroundRange(Cube *cube, Latitude &minLat,
      Latitude &maxLat, Longitude &minLon, Longitude &maxLon,
      bool allowEstimation) {
    // Do we need a RingRange method?
    // For now just return false
    if (HasCamera())
      if (p_camera->target()->shape()->name() == "Plane") return false;
    if (HasProjection())
      if (p_projection->projectionType() == Projection::RingPlane) return false;

    minLat = Latitude();
    maxLat = Latitude();
    minLon = Longitude();
    maxLon = Longitude();

    // If we have a footprint, use it
    try {
      if (cube) {
        ImagePolygon poly;
        cube->read(poly);
        geos::geom::MultiPolygon *footprint = PolygonTools::MakeMultiPolygon(
            poly.Polys()->clone());

        geos::geom::Geometry *envelope = footprint->getEnvelope();
        geos::geom::CoordinateSequence *coords = envelope->getCoordinates();

        for (unsigned int i = 0; i < coords->getSize(); i++) {
          const geos::geom::Coordinate &coord = coords->getAt(i);

          Latitude coordLat(coord.y, Angle::Degrees);
          Longitude coordLon(coord.x, Angle::Degrees);

          if (!minLat.isValid() || minLat > coordLat)
            minLat = coordLat;
          if (!maxLat.isValid() || maxLat < coordLat)
            maxLat = coordLat;

          if (!minLon.isValid() || minLon > coordLon)
            minLon = coordLon;
          if (!maxLon.isValid() || maxLon < coordLon)
            maxLon = coordLon;
        }

        delete coords;
        coords = NULL;

        delete envelope;
        envelope = NULL;

        delete footprint;
        footprint = NULL;
      }
    }
    catch (IException &) {
    }

    if (!minLat.isValid() || !maxLat.isValid() ||
        !minLon.isValid() || !maxLon.isValid()) {
      if (HasCamera()) {
        // Footprint failed, ask the camera
        PvlGroup mappingGrp("Mapping");
        mappingGrp += PvlKeyword("LatitudeType", "Planetocentric");
        mappingGrp += PvlKeyword("LongitudeDomain", "360");
        mappingGrp += PvlKeyword("LongitudeDirection", "PositiveEast");

        Pvl mappingPvl;
        mappingPvl += mappingGrp;
        double minLatDouble;
        double maxLatDouble;
        double minLonDouble;
        double maxLonDouble;
        p_camera->GroundRange(
            minLatDouble, maxLatDouble,
            minLonDouble, maxLonDouble, mappingPvl);
        minLat = Latitude(minLatDouble, Angle::Degrees);
        maxLat = Latitude(maxLatDouble, Angle::Degrees);
        minLon = Longitude(minLonDouble, Angle::Degrees);
        maxLon = Longitude(maxLonDouble, Angle::Degrees);
      }
      else if (HasProjection()) {
        // Footprint failed, look in the mapping group
        PvlGroup mappingGrp = p_projection->Mapping();
        if (mappingGrp.hasKeyword("MinimumLatitude") &&
            mappingGrp.hasKeyword("MaximumLatitude") &&
            mappingGrp.hasKeyword("MinimumLongitude") &&
            mappingGrp.hasKeyword("MaximumLongitude")) {

          minLat = Latitude(mappingGrp["MinimumLatitude"],
                            mappingGrp, Angle::Degrees);
          maxLat = Latitude(mappingGrp["MaximumLatitude"],
                            mappingGrp, Angle::Degrees);
          minLon = Longitude(mappingGrp["MinimumLongitude"],
                             mappingGrp, Angle::Degrees);
          maxLon = Longitude(mappingGrp["MaximumLongitude"],
                             mappingGrp, Angle::Degrees);

        }
        else if (allowEstimation && cube) {
          // Footprint and mapping failed... no lat/lon range of any kind is
          //   available. Let's test points in the image to try to make our own
          //   extent.
          QList<QPointF> imagePoints;

          // Reset to TProjection
          TProjection *tproj = (TProjection *) p_projection;

                    /*
           * This is where we're testing:
           *
           *  |---------------|
           *  |***************|
           *  |**     *     **|
           *  |*  *   *   *  *|
           *  |*    * * *    *|
           *  |***************|
           *  |*    * * *    *|
           *  |*  *   *   *  *|
           *  |**     *     **|
           *  |***************|
           *  |---------------|
           *
           * We'll test at the edges, a plus (+) and an (X) to help DEMs work.
           */

          int sampleCount = cube->sampleCount();
          int lineCount = cube->lineCount();

          int stepsPerLength = 20; //number of steps per length
          double aspectRatio = (double)lineCount / (double)sampleCount;
          double xStepSize = sampleCount / stepsPerLength;
          double yStepSize = xStepSize * aspectRatio;

          if (lineCount > sampleCount) {
            aspectRatio = (double)sampleCount / (double)lineCount;
            yStepSize = lineCount / stepsPerLength;
            xStepSize = yStepSize * aspectRatio;
          }

          double yWalked = 0.5;

          //3 vertical lines
          for (int i = 0; i < 3; i++) {
            double xValue = 0.5 + ( i * (sampleCount / 2) );

            while (yWalked <= lineCount) {
              imagePoints.append( QPointF(xValue, yWalked) );
              yWalked += yStepSize;
            }

            yWalked = 0.5;
          }

          double xWalked = 0.5;

          //3 horizontal lines
          for (int i = 0; i < 3; i++) {
            double yValue = 0.5 + ( i * (lineCount / 2) );

            while (xWalked <= sampleCount) {
              imagePoints.append( QPointF(xWalked, yValue) );
              xWalked += xStepSize;
            }

            xWalked = 0.5;
          }

          double xDiagonalWalked = 0.5;
          double yDiagonalWalked = 0.5;
          xStepSize = sampleCount / stepsPerLength;
          yStepSize = lineCount / stepsPerLength;

          //Top-Down Diagonal
          while ( (xDiagonalWalked <= sampleCount) && (yDiagonalWalked <= lineCount) ) {
            imagePoints.append( QPointF(xDiagonalWalked, yDiagonalWalked) );
            xDiagonalWalked += xStepSize;
            yDiagonalWalked += yStepSize;
          }

          xDiagonalWalked = 0.5;

          //Bottom-Up Diagonal
          while ( (xDiagonalWalked <= sampleCount) && (yDiagonalWalked >= 0) ) {
            imagePoints.append( QPointF(xDiagonalWalked, yDiagonalWalked) );
            xDiagonalWalked += xStepSize;
            yDiagonalWalked -= yStepSize;
          }

          foreach (QPointF imagePoint, imagePoints) {
            if (tproj->SetWorld(imagePoint.x(), imagePoint.y())) {
              Latitude latResult(tproj->UniversalLatitude(),
                                 Angle::Degrees);
              Longitude lonResult(tproj->UniversalLongitude(),
                                  Angle::Degrees);
              if (minLat.isValid())
                minLat = qMin(minLat, latResult);
              else
                minLat = latResult;

              if (maxLat.isValid())
                maxLat = qMax(maxLat, latResult);
              else
                maxLat = latResult;

              if (minLon.isValid())
                minLon = qMin(minLon, lonResult);
              else
                minLon = lonResult;

              if (maxLon.isValid())
                maxLon = qMax(maxLon, lonResult);
              else
                maxLon = lonResult;
            }
          }
        }
      }
    }

    return (minLat.isValid() && maxLat.isValid() &&
            minLon.isValid() && maxLon.isValid() &&
            minLat < maxLat && minLon < maxLon);
  }
}