Camera.cpp

Go to the documentation of this file.

#include "Camera.h"

#include <cfloat>
#include <cmath>
#include <iomanip>

#include <QDebug>
#include <QList>
#include <QPair>
#include <QTime>
#include <QVector>

#include "Angle.h"
#include "Constants.h"
#include "CameraDetectorMap.h"
#include "CameraFocalPlaneMap.h"
#include "CameraDistortionMap.h"
#include "CameraGroundMap.h"
#include "CameraSkyMap.h"
#include "DemShape.h"
#include "IException.h"
#include "IString.h"
#include "iTime.h"
#include "Latitude.h"
#include "Longitude.h"
#include "NaifStatus.h"
#include "Projection.h"
#include "ProjectionFactory.h"
#include "RingPlaneProjection.h"
#include "ShapeModel.h"
#include "SurfacePoint.h"
#include "Target.h"
#include "TProjection.h"

using namespace std;

namespace Isis {
  Camera::Camera(Cube &cube) : Sensor(cube) {
    // Get the image size which can be different than the alpha cube size
    p_lines = cube.lineCount();
    p_samples = cube.sampleCount();
    p_bands = cube.bandCount();

    SetGeometricTilingHint();

    // Get the AlphaCube information
    p_alphaCube = new AlphaCube(cube);

    // Get the projection group if it exists
    Pvl &lab = *cube.label();
    if (lab.findObject("IsisCube").hasGroup("Mapping")) {
      p_projection = ProjectionFactory::CreateFromCube(lab);
    }
    else {
      p_projection = NULL;
    }
    p_ignoreProjection = false;

    // Initialize stuff
    p_focalLength = 0.0;
    p_pixelPitch = 1.0;
    p_referenceBand = 0;
    p_childBand = 1;

    p_distortionMap = NULL;
    p_focalPlaneMap = NULL;
    p_detectorMap = NULL;
    p_groundMap = NULL;
    p_skyMap = NULL;

    // See if we have a reference band
    PvlGroup &inst = lab.findObject("IsisCube").findGroup("Instrument");
    if (inst.hasKeyword("ReferenceBand")) {
      p_referenceBand = inst["ReferenceBand"];
    }

    p_groundRangeComputed = false;
    p_raDecRangeComputed = false;
    p_ringRangeComputed = false;
    p_pointComputed = false;
  }

  Camera::~Camera() {
    if (p_projection) {
      delete p_projection;
      p_projection = NULL;
    }

    if (p_alphaCube) {
      delete p_alphaCube;
      p_alphaCube = NULL;
    }

    if (p_distortionMap) {
      delete p_distortionMap;
      p_distortionMap = NULL;
    }

    if (p_focalPlaneMap) {
      delete p_focalPlaneMap;
      p_focalPlaneMap = NULL;
    }

    if (p_detectorMap) {
      delete p_detectorMap;
      p_detectorMap = NULL;
    }

    if (p_groundMap) {
      delete p_groundMap;
      p_groundMap = NULL;
    }

    if (p_skyMap) {
      delete p_skyMap;
      p_skyMap = NULL;
    }
  }

  bool Camera::SetImage(const double sample, const double line) {
    p_childSample = sample;
    p_childLine = line;
    p_pointComputed = true;

    // get shape
    // TODO: we need to validate this pointer (somewhere)
    ShapeModel *shape = target()->shape();

    // Case of no map projection
    if (p_projection == NULL || p_ignoreProjection) {
      // Convert to parent coordinate (remove crop, pad, shrink, enlarge)
      double parentSample = p_alphaCube->AlphaSample(sample);
      double parentLine = p_alphaCube->AlphaLine(line);
      //cout << "\n\n\n\n\n\nCube: " << sample << " " << line << endl;
      //cout << "alpha cube: " << parentSample << " " << parentLine << endl; //debug
      //cout.precision(15);
      //cout << "Time: " << Time().Et() << endl;
      bool success = false;
      success = p_detectorMap->SetParent(parentSample, parentLine);
      // Convert from parent to detector
      if (success) {
        double detectorSample = p_detectorMap->DetectorSample();
        double detectorLine   = p_detectorMap->DetectorLine();
        // cout << "detector: " << detectorSample << " " << detectorLine << endl;  //debug
        success = p_focalPlaneMap->SetDetector(detectorSample, detectorLine);
        // Now Convert from detector to distorted focal plane
        if (success) {
          double focalPlaneX = p_focalPlaneMap->FocalPlaneX();
          double focalPlaneY = p_focalPlaneMap->FocalPlaneY();
          // cout << "focal plane: " << focalPlaneX << " " << focalPlaneY << endl; //debug
          success = p_distortionMap->SetFocalPlane(focalPlaneX, focalPlaneY);
          // Remove optical distortion
          if (success) {
            // Map to the ground
            double x = p_distortionMap->UndistortedFocalPlaneX();
            double y = p_distortionMap->UndistortedFocalPlaneY();
            double z = p_distortionMap->UndistortedFocalPlaneZ();
            return p_groundMap->SetFocalPlane(x, y, z);
          }
        }
      }
    }

    // The projection is a sky map
    else if (p_projection->IsSky()) {
      TProjection *tproj = (TProjection *) p_projection;
      if (tproj->SetWorld(sample, line)) {
        if (SetRightAscensionDeclination(tproj->Longitude(),
                                        tproj->UniversalLatitude())) {
          p_childSample = sample;
          p_childLine = line;

          return HasSurfaceIntersection();
        }
      }
    }

    // We have map projected camera model
    else {
      Latitude lat;
      Longitude lon;
      Distance rad;
      if (shape->name() != "Plane") { // this is the normal behavior
        if (p_projection->SetWorld(sample, line)) {
          TProjection *tproj = (TProjection *) p_projection;
          lat = Latitude(tproj->UniversalLatitude(), Angle::Degrees);
          lon = Longitude(tproj->UniversalLongitude(), Angle::Degrees);
          rad = Distance(LocalRadius(lat, lon));
          if (!rad.isValid()) {
            shape->setHasIntersection(false);
            return false;
          }
          SurfacePoint surfPt(lat, lon, rad);
          if (SetGround(surfPt)) {
            p_childSample = sample;
            p_childLine = line;

            shape->setHasIntersection(true);
            return true;
          }
        }
      }
      else { // shape is ring plane
        if (p_projection->SetWorld(sample, line)) {
          RingPlaneProjection *rproj = (RingPlaneProjection *) p_projection;
          lat = Latitude(0.0, Angle::Degrees);
          lon = Longitude(rproj->UniversalRingLongitude(), Angle::Degrees);
          rad = Distance(rproj->UniversalRingRadius(),Distance::Meters);

          if (!rad.isValid()) {
            shape->setHasIntersection(false);
            return false;
          }
          SurfacePoint surfPt(lat, lon, rad);
          if (SetGround(surfPt)) {
            p_childSample = sample;
            p_childLine = line;

            shape->setHasIntersection(true);
            return true;
          }
        }
      }
    }

    // failure
    shape->clearSurfacePoint();
    return false;
  }

  bool Camera::SetUniversalGround(const double latitude, const double longitude) {
    // Convert lat/lon or rad/az (i.e. ring rad / ring lon) to undistorted focal plane x/y
    if (p_groundMap->SetGround(Latitude(latitude, Angle::Degrees),
                              Longitude(longitude, Angle::Degrees))) {
      return RawFocalPlanetoImage();
    }

    target()->shape()->clearSurfacePoint();
    return false;
  }


  bool Camera::SetGround(Latitude latitude, Longitude longitude) {
    ShapeModel *shape = target()->shape();
    Distance localRadius;

    if (shape->name() != "Plane") { // this is the normal behavior
      localRadius = LocalRadius(latitude, longitude);
    }
    else {
      localRadius = Distance(latitude.degrees(),Distance::Kilometers);
      latitude = Latitude(0.,Angle::Degrees);
    }

      if (!localRadius.isValid()) {
      target()->shape()->clearSurfacePoint();
      return false;
    }

    return SetGround(SurfacePoint(latitude, longitude, localRadius));
  }



  bool Camera::SetGround(const SurfacePoint & surfacePt) {
    ShapeModel *shape = target()->shape();
    if (!surfacePt.Valid()) {
      shape->clearSurfacePoint();
      return false;
    }

    // Convert lat/lon to undistorted focal plane x/y
    if (p_groundMap->SetGround(surfacePt)) {
      return RawFocalPlanetoImage();
    }

    shape->clearSurfacePoint();
    return false;
  }


  bool Camera::RawFocalPlanetoImage() {
    double ux = p_groundMap->FocalPlaneX();
    double uy = p_groundMap->FocalPlaneY();
    
    // get shape
    // TODO: we need to validate this pointer (somewhere)
    ShapeModel *shape = target()->shape();
   
    //cout << "undistorted focal plane: " << ux << " " << uy << endl; //debug
    //cout.precision(15);
    //cout << "Backward Time: " << Time().Et() << endl;
    // Convert undistorted x/y to distorted x/y
    bool success = p_distortionMap->SetUndistortedFocalPlane(ux, uy);
    if (success) {
      double focalPlaneX = p_distortionMap->FocalPlaneX();
      double focalPlaneY = p_distortionMap->FocalPlaneY();
      //cout << "focal plane: " << focalPlaneX << " " << focalPlaneY << endl; //debug
      // Convert distorted x/y to detector position
      success = p_focalPlaneMap->SetFocalPlane(focalPlaneX, focalPlaneY);
      if (success) {
        double detectorSample = p_focalPlaneMap->DetectorSample();
        double detectorLine = p_focalPlaneMap->DetectorLine();
        //cout << "detector: " << detectorSample << " " << detectorLine << endl;
        // Convert detector to parent position
        success = p_detectorMap->SetDetector(detectorSample, detectorLine);
        if (success) {
          double parentSample = p_detectorMap->ParentSample();
          double parentLine = p_detectorMap->ParentLine();
          //cout << "cube: " << parentSample << " " << parentLine << endl; //debug
          if (p_projection == NULL || p_ignoreProjection) {
            p_childSample = p_alphaCube->BetaSample(parentSample);
            p_childLine = p_alphaCube->BetaLine(parentLine);
            p_pointComputed = true;
            shape->setHasIntersection(true);
            return true;
          }
          else if (p_projection->IsSky()) {
            if (p_projection->SetGround(Declination(), RightAscension())) {
              p_childSample = p_projection->WorldX();
              p_childLine = p_projection->WorldY();
              p_pointComputed = true;
              shape->setHasIntersection(true);
              return true;
            }
          }
          else  if (p_projection->projectionType() == Projection::Triaxial) {
            if (p_projection->SetUniversalGround(UniversalLatitude(), UniversalLongitude())) {
              p_childSample = p_projection->WorldX();
              p_childLine = p_projection->WorldY();
              p_pointComputed = true;
              shape->setHasIntersection(true);
              return true;
            }
          }
          else { // ring plane
            // UniversalLongitude should return azimuth (ring longitude) in this case TODO: when we make the change
            // to real azimuths this value may need to be adjusted or code changed in the shapemodel or
            // surfacepoint class.
            if (p_projection->SetUniversalGround(LocalRadius().meters(), UniversalLongitude())) {
              p_childSample = p_projection->WorldX();
              p_childLine = p_projection->WorldY();
              p_pointComputed = true;
              shape->setHasIntersection(true);
              return true;
            }
          }
        }
      }
    }

   shape->clearSurfacePoint();
   return false;
  }



  bool Camera::SetUniversalGround(const double latitude, const double longitude,
                                  const double radius) {
    // Convert lat/lon to undistorted focal plane x/y
    if (p_groundMap->SetGround(SurfacePoint(Latitude(latitude, Angle::Degrees),
                                            Longitude(longitude, Angle::Degrees),
                                            Distance(radius, Distance::Meters)))) {
      return RawFocalPlanetoImage();  // sets p_hasIntersection
    }

    target()->shape()->clearSurfacePoint();
   return false;
  }

  double Camera::DetectorResolution() {
    if (HasSurfaceIntersection()) {
      double sB[3];
      instrumentPosition(sB);
      double pB[3];
      Coordinate(pB);
      double a = sB[0] - pB[0];
      double b = sB[1] - pB[1];
      double c = sB[2] - pB[2];
      double dist = sqrt(a * a + b * b + c * c) * 1000.0;
      return dist / (p_focalLength / p_pixelPitch);
    }
    return -1.0;
  }

  double Camera::SampleResolution() {
    return DetectorResolution() * p_detectorMap->SampleScaleFactor();
  }

  double Camera::LineResolution() {
    return DetectorResolution() * p_detectorMap->LineScaleFactor();
  }

  double Camera::PixelResolution() {
    double lineRes = LineResolution();
    double sampRes = SampleResolution();
    if (lineRes < 0.0) return -1.0;
    if (sampRes < 0.0) return -1.0;
    return (lineRes + sampRes) / 2.0;
  }

  double Camera::LowestImageResolution() {
    GroundRangeResolution();
    return p_maxres;
  }

  double Camera::HighestImageResolution() {
    GroundRangeResolution();
    return p_minres;
  }

  void Camera::GroundRangeResolution() {
    // Software adjustment is needed if we get here -- call RingRangeResolution instead
    if (target()->shape()->name() == "Plane") {
      IString msg = "Images with plane targets should use Camera method RingRangeResolution ";
      msg += "instead of GroundRangeResolution";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // Have we already done this
    if (p_groundRangeComputed) return;
    p_groundRangeComputed = true;

    bool computed = p_pointComputed;
    double originalSample = Sample();
    double originalLine = Line();
    int originalBand = Band();

    // Initializations
    p_minlat    = DBL_MAX;
    p_minlon    = DBL_MAX;
    p_minlon180 = DBL_MAX;
    p_maxlat    = -DBL_MAX;
    p_maxlon    = -DBL_MAX;
    p_maxlon180 = -DBL_MAX;
    p_minres    = DBL_MAX;
    p_maxres    = -DBL_MAX;

    // See if we have band dependence and loop for the appropriate number of bands
    int eband = p_bands;
    if (IsBandIndependent()) eband = 1;
    for (int band = 1; band <= eband; band++) {
      SetBand(band);

      // Loop for each line testing the left and right sides of the image
      for (int line = 1; line <= p_lines + 1; line++) {
        // Look for the first good lat/lon on the left edge of the image
        // If it is the first or last line then test the whole line
        int samp;
        for (samp = 1; samp <= p_samples + 1; samp++) {

          if (SetImage((double)samp - 0.5, (double)line - 0.5)) {
            double lat = UniversalLatitude();
            double lon = UniversalLongitude();
            if (lat < p_minlat) p_minlat = lat;
            if (lat > p_maxlat) p_maxlat = lat;
            if (lon < p_minlon) p_minlon = lon;
            if (lon > p_maxlon) p_maxlon = lon;

            if (lon > 180.0) lon -= 360.0;
            if (lon < p_minlon180) p_minlon180 = lon;
            if (lon > p_maxlon180) p_maxlon180 = lon;

            double res = PixelResolution();
            if (res > 0.0) {
              if (res < p_minres) p_minres = res;
              if (res > p_maxres) p_maxres = res;
            }
            if ((line != 1) && (line != p_lines + 1)) break;
          }
        } // end loop through samples

        //We've already checked the first and last lines.
        if (line == 1) continue;
        if (line == p_lines + 1) continue;

        // Look for the first good lat/lon on the right edge of the image
        if (samp < p_samples + 1) {
          for (samp = p_samples + 1; samp >= 1; samp--) {
            if (SetImage((double)samp - 0.5, (double)line - 0.5)) {
              double lat = UniversalLatitude();
              double lon = UniversalLongitude();
              if (lat < p_minlat) p_minlat = lat;
              if (lat > p_maxlat) p_maxlat = lat;
              if (lon < p_minlon) p_minlon = lon;
              if (lon > p_maxlon) p_maxlon = lon;

              if (lon > 180.0) lon -= 360.0;
              if (lon < p_minlon180) p_minlon180 = lon;
              if (lon > p_maxlon180) p_maxlon180 = lon;

              double res = PixelResolution();
              if (res > 0.0) {
                if (res < p_minres) p_minres = res;
                if (res > p_maxres) p_maxres = res;
              }
              break;
            }
          }
        }
      } // end loop through lines

      // Test at the sub-spacecraft point to see if we have a
      // better resolution
      double lat, lon;

      subSpacecraftPoint(lat, lon);
      Latitude latitude(lat, Angle::Degrees);
      Longitude longitude(lon, Angle::Degrees);
      // get the local radius for the subspacecraft point
      Distance radius(LocalRadius(latitude, longitude));
      SurfacePoint testPoint;

      if (radius.isValid()) {

        testPoint = SurfacePoint(latitude, longitude, radius);

        if (SetGround(testPoint)) {
          if (Sample() >= 0.5 && Line() >= 0.5 &&
              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
            double res = PixelResolution();
            if (res > 0.0) {
              if (res < p_minres) p_minres = res;
              if (res > p_maxres) p_maxres = res;
            }
          }
        }
      } // end valid local (subspacecraft) radius

      // Special test for ground range to see if either pole is in the image
      latitude = Latitude(90, Angle::Degrees);
      longitude = Longitude(0.0, Angle::Degrees);
      // get radius for north pole
      radius = LocalRadius(latitude, longitude);

      if (radius.isValid()) {

        testPoint = SurfacePoint(latitude, longitude, radius);

        if (SetGround(testPoint)) {
          if (Sample() >= 0.5 && Line() >= 0.5 &&
              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
            p_maxlat = 90.0;
            p_minlon = 0.0;
            p_maxlon = 360.0;
            p_minlon180 = -180.0;
            p_maxlon180 = 180.0;
          }
        }
      } // end valid north polar radius

      latitude = Latitude(-90, Angle::Degrees);
      // get radius for south pole
      radius = LocalRadius(latitude, longitude);

      if (radius.isValid()) {

        testPoint = SurfacePoint(latitude, longitude, radius);
        if (SetGround(testPoint)) {
          if (Sample() >= 0.5 && Line() >= 0.5 &&
              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
            p_minlat = -90.0;
            p_minlon = 0.0;
            p_maxlon = 360.0;
            p_minlon180 = -180.0;
            p_maxlon180 = 180.0;
          }
        }
      } // end valid south polar radius

      // Another special test for ground range as we could have the
      // 0-360 seam running right through the image so
      // test it as well (the increment may not be fine enough !!!)
      for (Latitude lat = Latitude(p_minlat, Angle::Degrees);
                    lat <= Latitude(p_maxlat, Angle::Degrees);
                    lat += Angle((p_maxlat - p_minlat) / 10.0, Angle::Degrees)) {
        if (SetGround(lat, Longitude(0.0, Angle::Degrees))) {
          if (Sample() >= 0.5 && Line() >= 0.5 &&
              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
            p_minlon = 0.0;
            p_maxlon = 360.0;
            break;
          }
        } // end if set ground (lat, 0)

        // Another special test for ground range as we could have the
        // -180-180 seam running right through the image so
        // test it as well (the increment may not be fine enough !!!)
        if (SetGround(lat, Longitude(180.0, Angle::Degrees))) {
          if (Sample() >= 0.5 && Line() >= 0.5 &&
              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
            p_minlon180 = -180.0;
            p_maxlon180 = 180.0;
            break;
          }
        } // end if set ground (lat, 180)
      } // end for loop (latitudes from min to max)
    } // end for loop through bands

    SetBand(originalBand);

    if(computed) {
      SetImage(originalSample, originalLine);
    }
    else {
      p_pointComputed = false;
    }

    if (p_minlon == DBL_MAX  ||  p_minlat == DBL_MAX  ||  p_maxlon == -DBL_MAX  ||  p_maxlat == -DBL_MAX) {
      string message = "Camera missed planet or SPICE data off.";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }


    // Checks for invalid lat/lon ranges
//    if(p_minlon == DBL_MAX  ||  p_minlat == DBL_MAX  ||  p_maxlon == -DBL_MAX  ||  p_maxlat == -DBL_MAX) {
//      string message = "Camera missed planet or SPICE data off.";
//      throw IException(IException::Unknown, message, _FILEINFO_);
//    }
  }


  void Camera::ringRangeResolution() {
    // TODO Add test to make sure we have a ring plane image **

    // Have we already done this
    if (p_ringRangeComputed) return;

    p_ringRangeComputed = true;

    bool computed = p_pointComputed;
    double originalSample = Sample();
    double originalLine = Line();
    int originalBand = Band();

    // Initializations
    p_minRingRadius       =  DBL_MAX;
    p_minRingLongitude    =  DBL_MAX;
    p_minRingLongitude180 =  DBL_MAX;
    p_maxRingRadius       = -DBL_MAX;
    p_maxRingLongitude    = -DBL_MAX;
    p_maxRingLongitude180 = -DBL_MAX;
    p_minres              =  DBL_MAX;
    p_maxres              = -DBL_MAX;

    // See if we have band dependence and loop for the appropriate number of bands
    int eband = p_bands;
    if (IsBandIndependent())
      eband = 1;

    for (int band = 1; band <= eband; band++) {
      SetBand(band);

      // Loop for each line testing the left and right sides of the image
      for (int line = 1; line <= p_lines + 1; line++) {

        // Look for the first good radius/azimuth on the left edge of the image
        // If it is the first or last line then test the whole line
        int samp;
        for (samp = 1; samp <= p_samples + 1; samp++) {

          if (SetImage((double)samp - 0.5, (double)line - 0.5)) {
            double radius = LocalRadius().meters();
            double azimuth = UniversalLongitude();
            if (radius < p_minRingRadius) p_minRingRadius = radius;
            if (radius > p_maxRingRadius) p_maxRingRadius = radius;
            if (azimuth < p_minRingLongitude) p_minRingLongitude = azimuth;
            if (azimuth > p_maxRingLongitude) p_maxRingLongitude = azimuth;

            if (azimuth > 180.0) azimuth -= 360.0;
            if (azimuth < p_minRingLongitude180) p_minRingLongitude180 = azimuth;
            if (azimuth > p_maxRingLongitude180) p_maxRingLongitude180 = azimuth;

            double res = PixelResolution();
            if (res > 0.0) {
              if (res < p_minres) p_minres = res;
              if (res > p_maxres) p_maxres = res;
            }
            if ((line != 1) && (line != p_lines + 1)) break;
          }
        }
        //We've already checked the first and last lines.
        if (line == 1) continue;
        if (line == p_lines + 1) continue;

        // Look for the first good rad/azimuth on the right edge of the image
        if (samp < p_samples + 1) {
          for(samp = p_samples + 1; samp >= 1; samp--) {
            if (SetImage((double)samp - 0.5, (double)line - 0.5)) {
              double radius = LocalRadius().meters();
              double azimuth = UniversalLongitude();
              if (radius < p_minRingRadius) p_minRingRadius = radius;
              if (radius > p_maxRingRadius) p_maxRingRadius = radius;
              if (azimuth < p_minRingLongitude) p_minRingLongitude = azimuth;
              if (azimuth > p_maxRingLongitude) p_maxRingLongitude = azimuth;

              if (azimuth > 180.0) azimuth -= 360.0;
              if (azimuth < p_minRingLongitude180) p_minRingLongitude180 = azimuth;
              if (azimuth > p_maxRingLongitude180) p_maxRingLongitude180 = azimuth;

              double res = PixelResolution();
              if (res > 0.0) {
                if (res < p_minres) p_minres = res;
                if (res > p_maxres) p_maxres = res;
              }
              break;
            }
          }
        }
      }

      // Test at the sub-spacecraft point to see if we have a
      // better resolution
      // TODO: is there something to this analogous for ring images?
//      double rad, lon;

//      subSpacecraftPoint(lat, lon);
//      Latitude latitude(lat, Angle::Degrees);
//      Longitude longitude(lon, Angle::Degrees);
//      Distance radius(LocalRadius(latitude, longitude));
//      SurfacePoint testPoint;

//      if (radius.isValid()) {

//        testPoint = SurfacePoint(latitude, longitude, radius);

//        if(SetGround(testPoint)) {
//          if(Sample() >= 0.5 && Line() >= 0.5 &&
//              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
//            double res = PixelResolution();
//            if(res > 0.0) {
//              if(res < p_minres) p_minres = res;
//              if(res > p_maxres) p_maxres = res;
//            }
//          }
//        }
//      }

      // Another special test for ring range as we could have the
      // 0-360 seam running right through the image so
      // test it as well (the increment may not be fine enough !!!)
      for (Distance radius = Distance(p_minRingRadius, Distance::Meters);
                   radius <= Distance(p_maxRingRadius, Distance::Meters);
                   radius += Distance((p_maxRingRadius - p_minRingRadius) / 10.0, Distance::Meters)) {
        if (SetGround(SurfacePoint(Latitude(0.0, Angle::Degrees), Longitude(0.0, Angle::Degrees), radius))) {
          if (Sample() >= 0.5 && Line() >= 0.5 &&
              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
            p_minRingLongitude = 0.0;
            p_maxRingLongitude = 360.0;
            break;
          }
        }

      // for (Latitude lat = Latitude(p_minlat, Angle::Degrees);
      //              lat <= Latitude(p_maxlat, Angle::Degrees);
      //              lat += Angle((p_maxlat - p_minlat) / 10.0, Angle::Degrees)) {
      //   if (SetGround(lat, Longitude(0.0, Angle::Degrees))) {
      //     if (Sample() >= 0.5 && Line() >= 0.5 &&
      //         Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
      //       p_minlon = 0.0;
      //       p_maxlon = 360.0;
      //       break;
      //     }
      //   }

        // Another special test for ring range as we could have the
        // -180-180 seam running right through the image so
        // test it as well (the increment may not be fine enough !!!)
        if (SetGround(Latitude(0.0, Angle::Degrees), Longitude(180.0, Angle::Degrees))) {
          if (Sample() >= 0.5 && Line() >= 0.5 &&
              Sample() <= p_samples + 0.5 && Line() <= p_lines + 0.5) {
            p_minRingLongitude180 = -180.0;
            p_maxRingLongitude180 = 180.0;
            break;
          }
        }
      }
    } // end loop over bands

    SetBand(originalBand);

    if (computed) {
      SetImage(originalSample, originalLine);
    }
    else {
      p_pointComputed = false;
    }

    // Checks for invalid radius/lon ranges
    if (p_minRingRadius == DBL_MAX  ||  p_minRingRadius == DBL_MAX  || p_minRingLongitude == DBL_MAX  ||  p_maxRingLongitude == -DBL_MAX) {
      string message = "RingPlane ShapeModel - Camera missed plane or SPICE data off.";
      throw IException(IException::Unknown, message, _FILEINFO_);
    }
  }


  bool Camera::IntersectsLongitudeDomain(Pvl &pvl) {
    double minlat, minlon, maxlat, maxlon;
    return GroundRange(minlat, maxlat, minlon, maxlon, pvl);
  }

  bool Camera::GroundRange(double &minlat, double &maxlat,
                           double &minlon, double &maxlon,
                           Pvl &pvl) {
    // Compute the ground range and resolution
    GroundRangeResolution();

    // Get the default radii
    Distance localRadii[3];
    radii(localRadii);
    Distance &a = localRadii[0];
    Distance &b = localRadii[2];

    // See if the PVL overrides the radii
    PvlGroup map = pvl.findGroup("Mapping", Pvl::Traverse);

    if(map.hasKeyword("EquatorialRadius"))
      a = Distance(toDouble(map["EquatorialRadius"][0]), Distance::Meters);

    if(map.hasKeyword("PolarRadius"))
      b = Distance(toDouble(map["PolarRadius"][0]), Distance::Meters);

    // Convert to planetographic if necessary
    minlat = p_minlat;
    maxlat = p_maxlat;
    if(map.hasKeyword("LatitudeType")) {
      QString latType = (QString) map["LatitudeType"];
      if (latType.toUpper() == "PLANETOGRAPHIC") {
        if (abs(minlat) < 90.0) {  // So tan doesn't fail
          minlat *= PI / 180.0;
          minlat = atan(tan(minlat) * (a / b) * (a / b));
          minlat *= 180.0 / PI;
        }

        if(abs(maxlat) < 90.0) {  // So tan doesn't fail
          maxlat *= PI / 180.0;
          maxlat = atan(tan(maxlat) * (a / b) * (a / b));
          maxlat *= 180.0 / PI;
        }
      }
    }

    // Assume 0 to 360 domain but change it if necessary
    minlon = p_minlon;
    maxlon = p_maxlon;
    bool domain360 = true;
    if(map.hasKeyword("LongitudeDomain")) {
      QString lonDomain = (QString) map["LongitudeDomain"];
      if(lonDomain.toUpper() == "180") {
        minlon = p_minlon180;
        maxlon = p_maxlon180;
        domain360 = false;
      }
    }

    // Convert to the proper longitude direction
    if(map.hasKeyword("LongitudeDirection")) {
      QString lonDirection = (QString) map["LongitudeDirection"];
      if(lonDirection.toUpper() == "POSITIVEWEST") {
        double swap = minlon;
        minlon = -maxlon;
        maxlon = -swap;
      }
    }

    // Convert to the proper longitude domain
    if(domain360) {
      while(minlon < 0.0) {
        minlon += 360.0;
        maxlon += 360.0;
      }
      while(minlon > 360.0) {
        minlon -= 360.0;
        maxlon -= 360.0;
      }
    }
    else {
      while(minlon < -180.0) {
        minlon += 360.0;
        maxlon += 360.0;
      }
      while(minlon > 180.0) {
        minlon -= 360.0;
        maxlon -= 360.0;
      }
    }

    // Now return if it crosses the longitude domain boundary
    if((maxlon - minlon) > 359.0) return true;
    return false;
  }

  bool Camera::ringRange(double &minRingRadius, double &maxRingRadius, 
                         double &minRingLongitude, double &maxRingLongitude, Pvl &pvl) {
    // Compute the ring range and resolution
    ringRangeResolution();

    // Get the mapping group
    PvlGroup map = pvl.findGroup("Mapping", Pvl::Traverse);

    // Get the ring radius range
    minRingRadius = p_minRingRadius;
    maxRingRadius = p_maxRingRadius;

    // Assume 0 to 360 domain but change it if necessary
    minRingLongitude = p_minRingLongitude;
    maxRingLongitude = p_maxRingLongitude;
    bool domain360 = true;
    if (map.hasKeyword("RingLongitudeDomain")) {
      QString ringLongitudeDomain = (QString) map["RingLongitudeDomain"];
      if (ringLongitudeDomain == "180") {
        minRingLongitude = p_minRingLongitude180;
        maxRingLongitude = p_maxRingLongitude180;
        domain360 = false;
      }
    }

    // Convert to the proper azimuth direction
    if (map.hasKeyword("RingLongitudeDirection")) {
      QString ringLongitudeDirection = (QString) map["RingLongitudeDirection"];
      if (ringLongitudeDirection.toUpper() == "Clockwise") {
        double swap = minRingLongitude;
        minRingLongitude = -maxRingLongitude;
        maxRingLongitude = -swap;
      }
    }

    // Convert to the proper azimuth domain
    if (domain360) {
      while (minRingLongitude < 0.0) {
        minRingLongitude += 360.0;
        maxRingLongitude += 360.0;
      }
      while (minRingLongitude > 360.0) {
        minRingLongitude -= 360.0;
        maxRingLongitude -= 360.0;
      }
    }
    else {
      while (minRingLongitude < -180.0) {
        minRingLongitude += 360.0;
        maxRingLongitude += 360.0;
      }
      while (minRingLongitude > 180.0) {
        minRingLongitude -= 360.0;
        maxRingLongitude -= 360.0;
      }
    }

    // Now return if it crosses the azimuth domain boundary
    if ((maxRingLongitude - minRingLongitude) > 359.0) {
      return true;
    }
    return false;
  }


  void Camera::BasicMapping(Pvl &pvl) {
    PvlGroup map("Mapping");
    map += PvlKeyword("TargetName", target()->name());

    std::vector<Distance> radii = target()->radii();
    map += PvlKeyword("EquatorialRadius", toString(radii[0].meters()), "meters");
    map += PvlKeyword("PolarRadius", toString(radii[2].meters()), "meters");
    
    map += PvlKeyword("LatitudeType", "Planetocentric");
    map += PvlKeyword("LongitudeDirection", "PositiveEast");
    map += PvlKeyword("LongitudeDomain", "360");

    GroundRangeResolution();
    map += PvlKeyword("MinimumLatitude", toString(p_minlat));
    map += PvlKeyword("MaximumLatitude", toString(p_maxlat));
    map += PvlKeyword("MinimumLongitude", toString(p_minlon));
    map += PvlKeyword("MaximumLongitude", toString(p_maxlon));
    map += PvlKeyword("PixelResolution", toString(p_minres));

    map += PvlKeyword("ProjectionName", "Sinusoidal");
    pvl.addGroup(map);
  }


  void Camera::basicRingMapping(Pvl &pvl) {
    if (target()->shape()->name() != "Plane") { // If we get here and we don't have a plane, throw an error
      IString msg = "A ring plane projection has been requested on an image whose shape is not a ring plane.  ";
      msg += "Rerun spiceinit with shape=RINGPLANE.  ";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    PvlGroup map("Mapping");
    map += PvlKeyword("TargetName", target()->name());

    map += PvlKeyword("RingLongitudeDirection", "CounterClockwise");
    map += PvlKeyword("RingLongitudeDomain", "360");

    ringRangeResolution();
    map += PvlKeyword("MinimumRingRadius", toString(p_minRingRadius));
    map += PvlKeyword("MaximumRingRadius", toString(p_maxRingRadius));
    map += PvlKeyword("MinimumRingLongitude", toString(p_minRingLongitude));
    map += PvlKeyword("MaximumRingLongitude", toString(p_maxRingLongitude));
    map += PvlKeyword("PixelResolution", toString(p_minres));

    map += PvlKeyword("ProjectionName", "Planar");
    pvl.addGroup(map);
  }


  void Camera::SetFocalLength() {
    int code = naifIkCode();
    QString key = "INS" + toString(code) + "_FOCAL_LENGTH";
    SetFocalLength(Spice::getDouble(key));
  }

  void Camera::SetPixelPitch() {
    int code = naifIkCode();
    QString key = "INS" + toString(code) + "_PIXEL_PITCH";
    SetPixelPitch(Spice::getDouble(key));
  }

  bool Camera::SetRightAscensionDeclination(const double ra, const double dec) {
    if (p_skyMap->SetSky(ra, dec)) {
      double ux = p_skyMap->FocalPlaneX();
      double uy = p_skyMap->FocalPlaneY();
      if (p_distortionMap->SetUndistortedFocalPlane(ux, uy)) {
        double dx = p_distortionMap->FocalPlaneX();
        double dy = p_distortionMap->FocalPlaneY();
        if (p_focalPlaneMap->SetFocalPlane(dx, dy)) {
          double detectorSamp = p_focalPlaneMap->DetectorSample();
          double detectorLine = p_focalPlaneMap->DetectorLine();
          if (p_detectorMap->SetDetector(detectorSamp, detectorLine)) {
            double parentSample = p_detectorMap->ParentSample();
            double parentLine = p_detectorMap->ParentLine();

            if (p_projection == NULL || p_ignoreProjection) {
              p_childSample = p_alphaCube->BetaSample(parentSample);
              p_childLine = p_alphaCube->BetaLine(parentLine);
              p_pointComputed = true;
              return true;
            }
            else if (p_projection->IsSky()) {
              if (p_projection->SetGround(dec, ra)) {
                p_childSample = p_projection->WorldX();
                p_childLine = p_projection->WorldY();
                p_pointComputed = true;
                return true;
              }
            }
            else if (target()->shape()->hasIntersection()) {
              if (p_projection->SetUniversalGround(UniversalLatitude(),
                                                  UniversalLongitude())) {
                p_childSample = p_projection->WorldX();
                p_childLine = p_projection->WorldY();
                p_pointComputed = true;
                return true;
              }
            }
          }
        }
      }
    }

    return false;
  }


  void Camera::GetLocalNormal(double normal[3]) {

    ShapeModel *shapeModel = target()->shape();
    if ( !shapeModel->hasIntersection()) {
      // if the shape is not intersected, then clearly there is no normal
      normal[0] = normal[1] = normal[2] = 0.0;
      return;
    }

    QVector<double *> cornerNeighborPoints(4);
    bool computed = p_pointComputed;
    if (QString::compare(shapeModel->name(), "DemShape", Qt::CaseInsensitive) == 0) {
      // If the shape model is not DEM, we don't use the neighbor points.
      // So, to avoid the potetially expesive SetImage() calls used to find the neighbors,
      // just pass in empty vector to shape model's calculateLocalNormal() method
      shapeModel->calculateLocalNormal(cornerNeighborPoints);
    }
    else { // attempt to find local normal for DEM shapes using 4 surrounding points on the image

      // As documented in the doxygen above, the goal of this method is to
      // calculate a normal vector to the surface using the 4 corner surrounding points.
      double samp = Sample();
      double line = Line();

      // order of points in vector is top, bottom, left, right
      QList< QPair< double, double > > surroundingPoints;
      surroundingPoints.append(qMakePair(samp, line - 0.5));
      surroundingPoints.append(qMakePair(samp, line + 0.5));
      surroundingPoints.append(qMakePair(samp - 0.5, line));
      surroundingPoints.append(qMakePair(samp + 0.5, line));

      // save input state to be restored on return
      double originalSample = samp;
      double originalLine = line;

      // now we have all four points in the image, so find the same points on the surface
      for (int i = 0; i < cornerNeighborPoints.size(); i++) {
        cornerNeighborPoints[i] = new double[3];
      }

      Latitude lat;
      Longitude lon;
      Distance radius;

      // if this is a dsk, we only need to use the existing intercept point (plate) normal then return 
      for (int i = 0; i < cornerNeighborPoints.size(); i++) {
        // If a surrounding point fails, set it to the original point
        if (!(SetImage(surroundingPoints[i].first, surroundingPoints[i].second))) {
          surroundingPoints[i].first = samp;
          surroundingPoints[i].second = line;

          // If the original point fails too, we can't get a normal.  Clean up and return.
          if (!(SetImage(surroundingPoints[i].first, surroundingPoints[i].second))) {
            normal[0] = 0.;
            normal[1] = 0.;
            normal[2] = 0.;

            // restore input state
            if (computed) {
              SetImage(originalSample, originalLine);
            } 
            else {
              p_pointComputed = false;
            }

            // free memory
            for (int i = 0; i < cornerNeighborPoints.size(); i++) {
              delete [] cornerNeighborPoints[i];
            }

            return;
          }
        }

        SurfacePoint surfacePoint = GetSurfacePoint();
        lat = surfacePoint.GetLatitude();
        lon = surfacePoint.GetLongitude();
        radius = LocalRadius(lat, lon);

        latrec_c(radius.kilometers(), lon.radians(), lat.radians(), cornerNeighborPoints[i]);
      }

      // if the first 2 surrounding points match or the last 2 surrounding points match, 
      // we can't get a normal.  Clean up and return. 
      if ((surroundingPoints[0].first == surroundingPoints[1].first &&
          surroundingPoints[0].second == surroundingPoints[1].second) ||
         (surroundingPoints[2].first == surroundingPoints[3].first &&
          surroundingPoints[2].second == surroundingPoints[3].second)) {
        normal[0] = 0.;
        normal[1] = 0.;
        normal[2] = 0.;

        // restore input state
        if (!computed) {
          SetImage(originalSample, originalLine);
        } 
        else {
          p_pointComputed = false;
        }

        // free memory
        for (int i = 0; i < cornerNeighborPoints.size(); i++)
          delete [] cornerNeighborPoints[i];

        return;
      }

      // Restore input state to original point before calculating normal
      SetImage(originalSample, originalLine);
      shapeModel->calculateLocalNormal(cornerNeighborPoints);

      // free memory
      for (int i = 0; i < cornerNeighborPoints.size(); i++)
        delete [] cornerNeighborPoints[i];
   
    }

    // restore input state if calculation failed and clean up.
    if (!shapeModel->hasNormal()) {
      p_pointComputed = false;
      return;
    }

    // restore failed computed state
    if (!computed) {
      p_pointComputed = false;
    }

    // Set the method normal values
    std::vector<double> localNormal(3);
    localNormal = shapeModel->normal();
    memcpy(normal, (double *) &localNormal[0], sizeof(double) * 3);
  }


  void Camera::LocalPhotometricAngles(Angle & phase, Angle & incidence,
      Angle & emission, bool &success) {

    // get local normal vector
    double normal[3];
    GetLocalNormal(normal);
    success = true;

    // Check to make sure normal is valid
    SpiceDouble mag;
    unorm_c(normal,normal,&mag);
    if (mag == 0.) {
      success = false;
      return;
    }

    // get a normalized surface spacecraft vector
    SpiceDouble surfSpaceVect[3], unitizedSurfSpaceVect[3], dist;
    std::vector<double> sB = bodyRotation()->ReferenceVector(
        instrumentPosition()->Coordinate());

    SpiceDouble pB[3];
    SurfacePoint surfacePoint = GetSurfacePoint();
    pB[0] = surfacePoint.GetX().kilometers();
    pB[1] = surfacePoint.GetY().kilometers();
    pB[2] = surfacePoint.GetZ().kilometers();

    vsub_c((SpiceDouble *) &sB[0], pB, surfSpaceVect);
    unorm_c(surfSpaceVect, unitizedSurfSpaceVect, &dist);

    // get a normalized surface sun vector
    SpiceDouble surfaceSunVect[3];
    vsub_c(m_uB, pB, surfaceSunVect);
    SpiceDouble unitizedSurfSunVect[3];
    unorm_c(surfaceSunVect, unitizedSurfSunVect, &dist);

    // use normalized surface spacecraft and surface sun vectors to calculate
    // the phase angle (in radians)
    phase = Angle(vsep_c(unitizedSurfSpaceVect, unitizedSurfSunVect),
        Angle::Radians);

    // use normalized surface spacecraft and local normal vectors to calculate
    // the emission angle (in radians)
    emission = Angle(vsep_c(unitizedSurfSpaceVect, normal),
        Angle::Radians);

    // use normalized surface sun and normal vectors to calculate the incidence
    // angle (in radians)
    incidence = Angle(vsep_c(unitizedSurfSunVect, normal),
        Angle::Radians);
  }


  bool Camera::RaDecRange(double &minra, double &maxra,
                          double &mindec, double &maxdec) {
    TProjection *tproj = (TProjection *) p_projection;
    if (p_projection != NULL && !tproj->IsSky()) {
      IString msg = "Camera::RaDecRange can not calculate a right ascension, declination range";
      msg += " for projected images which are not projected to sky";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    bool computed = p_pointComputed;
    double originalSample = Sample();
    double originalLine = Line();
    int originalBand = Band();

    // Have we already done this
    if (!p_raDecRangeComputed) {
      p_raDecRangeComputed = true;

      // Initializations
      p_mindec    = DBL_MAX;
      p_minra     = DBL_MAX;
      p_minra180  = DBL_MAX;
      p_maxdec    = -DBL_MAX;
      p_maxra     = -DBL_MAX;
      p_maxra180  = -DBL_MAX;

      // See if we have band dependence and loop for the appropriate number of bands
      int eband = p_bands;
      if (IsBandIndependent()) eband = 1;
      for (int band = 1; band <= eband; band++) {
        this->SetBand(band);

        for (int line = 1; line <= p_lines; line++) {
          // Test left, top, and bottom sides
          int samp;
          for (samp = 1; samp <= p_samples; samp++) {
            SetImage((double)samp, (double)line);
            double ra = RightAscension();
            double dec = Declination();
            if (ra < p_minra) p_minra = ra;
            if (ra > p_maxra) p_maxra = ra;
            if (dec < p_mindec) p_mindec = dec;
            if (dec > p_maxdec) p_maxdec = dec;

            if (ra > 180.0) ra -= 360.0;
            if (ra < p_minra180) p_minra180 = ra;
            if (ra > p_maxra180) p_maxra180 = ra;

            if ((line != 1) && (line != p_lines)) break;
          }

          // Test right side
          if (samp < p_samples) {
            for (samp = p_samples; samp >= 1; samp--) {
              SetImage((double)samp, (double)line);
              double ra = RightAscension();
              double dec = Declination();
              if (ra < p_minra) p_minra = ra;
              if (ra > p_maxra) p_maxra = ra;
              if (dec < p_mindec) p_mindec = dec;
              if (dec > p_maxdec) p_maxdec = dec;

              if (ra > 180.0) ra -= 360.0;
              if (ra < p_minra180) p_minra180 = ra;
              if (ra > p_maxra180) p_maxra180 = ra;

              break;
            }
          }
        }

        // Special test for ground range to see if either pole is in the image
        if (SetRightAscensionDeclination(0.0, 90.0)) {
          if ((Line() >= 0.5) && (Line() <= p_lines) &&
              (Sample() >= 0.5) && (Sample() <= p_samples)) {
            p_maxdec = 90.0;
            p_minra = 0.0;
            p_maxra = 360.0;
            p_minra180 = -180.0;
            p_maxra180 = 180.0;
          }
        }

        if (SetRightAscensionDeclination(0.0, -90.0)) {
          if ((Line() >= 0.5) && (Line() <= p_lines) &&
              (Sample() >= 0.5) && (Sample() <= p_samples)) {
            p_mindec = -90.0;
            p_minra = 0.0;
            p_maxra = 360.0;
            p_minra180 = -180.0;
            p_maxra180 = 180.0;
          }
        }

        // Another special test for ground range as we could have the
        // 0-360 seam running right through the image so
        // test it as well (the increment may not be fine enough !!!)
        for (double dec = p_mindec; dec <= p_maxdec; dec += (p_maxdec - p_mindec) / 10.0) {
          if (SetRightAscensionDeclination(0.0, dec)) {
            if ((Line() >= 0.5) && (Line() <= p_lines) &&
                (Sample() >= 0.5) && (Sample() <= p_samples)) {
              p_minra = 0.0;
              p_maxra = 360.0;
              break;
            }
          }
        }

        // Another special test for ground range as we could have the
        // 0-360 seam running right through the image so
        // test it as well (the increment may not be fine enough !!!)
        for (double dec = p_mindec; dec <= p_maxdec; dec += (p_maxdec - p_mindec) / 10.0) {
          if (SetRightAscensionDeclination(180.0, dec)) {
            if ((Line() >= 0.5) && (Line() <= p_lines) &&
                (Sample() >= 0.5) && (Sample() <= p_samples)) {
              p_minra180 = -180.0;
              p_maxra180 = 180.0;
              break;
            }
          }
        }
      }
    }

    minra = p_minra;
    maxra = p_maxra;
    mindec = p_mindec;
    maxdec = p_maxdec;

    SetBand(originalBand);

    if (computed) {
      SetImage(originalSample, originalLine);
    }
    else {
      p_pointComputed = false;
    }

    return true;
  }


  double Camera::RaDecResolution() {
    TProjection *tproj = (TProjection *) p_projection;
    if (p_projection != NULL && !tproj->IsSky()) {
      IString msg = "Camera::RaDecResolution can not calculate a right ascension, declination resolution";
      msg += " for projected images which are not projected to sky";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    bool computed = p_pointComputed;
    double originalSample = Sample();
    double originalLine = Line();
    int originalBand = Band();

    SetImage(1.0, 1.0);
    double ra1 = RightAscension();
    double dec1 = Declination();

    SetImage(1.0, (double)p_lines);
    double ra2 = RightAscension();
    double dec2 = Declination();

    double dist = (ra1 - ra2) * (ra1 - ra2) + (dec1 - dec2) * (dec1 - dec2);
    dist = sqrt(dist);
    double lineRes = dist / (p_lines - 1);

    SetImage((double)p_samples, 1.0);
    ra2 = RightAscension();
    dec2 = Declination();

    dist = (ra1 - ra2) * (ra1 - ra2) + (dec1 - dec2) * (dec1 - dec2);
    dist = sqrt(dist);
    double sampRes = dist / (p_samples - 1);

    SetBand(originalBand);

    if (computed) {
      SetImage(originalSample, originalLine);
    }
    else {
      p_pointComputed = false;
    }

    return (sampRes < lineRes) ? sampRes : lineRes;
  }

  double Camera::NorthAzimuth() {
    if (target()->shape()->name() == "Plane") {
      QString msg = "North Azimuth is not available for plane target shapes.";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    // Get the latitude of your current location using the shape model
    // specified in the image Kernels
    double lat = UniversalLatitude();
    // We are in northern hemisphere
    if (lat >= 0.0) {
      return ComputeAzimuth(LocalRadius(90.0, 0.0), 90.0, 0.0);
    }
    // We are in southern hemisphere
    else {
      double azimuth = ComputeAzimuth(LocalRadius(-90.0, 0.0),
                                      -90.0, 0.0) + 180.0;
      if (azimuth > 360.0) azimuth = azimuth - 360.0;
      return azimuth;
    }
  }

  double Camera::SunAzimuth() {
    double lat, lon;
    subSolarPoint(lat, lon);
    return ComputeAzimuth(LocalRadius(lat, lon), lat, lon);
  }

  double Camera::SpacecraftAzimuth() {
    double lat, lon;
    subSpacecraftPoint(lat, lon);
    return ComputeAzimuth(LocalRadius(lat, lon), lat, lon);
  }

  double Camera::ComputeAzimuth(Distance radius,
                                const double lat, const double lon) {
    // Make sure we are on the planet, if not, north azimuth is meaningless
    if (!HasSurfaceIntersection()) return Isis::Null;

    // Need to save the "state" of the camera so we can restore it when the
    // method is done
    bool computed = p_pointComputed;

    NaifStatus::CheckErrors();

    // Get the azimuth's origin point (the current position) and its radius
    SpiceDouble azimuthOrigin[3];
    Coordinate(azimuthOrigin);
    Distance originRadius = LocalRadius();
    if (!originRadius.isValid()) {
      return Isis::Null;
    }

    // Convert the point of interest to rectangular coordinates (x,y,z) in the
    // body-fixed coordinate system and use the azimuth's origin radius to avoid the
    // situation where the DEM does not cover the entire planet
    SpiceDouble pointOfInterestFromBodyCenter[3];
    latrec_c(originRadius.kilometers(), lon * PI / 180.0,
             lat * PI / 180.0, pointOfInterestFromBodyCenter);

    // Get the difference vector with its tail at the azimuth origin and its
    // head at the point of interest by subtracting vectors the vectors that
    // originate at the body center
    // 
    // pointOfInterest = pointOfInterestFromBodyCenter - azimuthOriginFromBodyCenter
    // 
    SpiceDouble pointOfInterest[3];
    vsub_c(pointOfInterestFromBodyCenter, azimuthOrigin, pointOfInterest);

    // Get the component of the difference vector pointOfInterestFromAzimuthOrigin that is
    // perpendicular to the origin point (i.e. project perpendicularly onto the reference plane).
    // This will result in a point of interest vector that is in the plane tangent to the surface at
    // the origin point
    SpiceDouble pointOfInterestProj[3];
    vperp_c(pointOfInterest, azimuthOrigin, pointOfInterestProj);

    // Unitize the tangent vector to a 1 km length vector
    SpiceDouble pointOfInterestProjUnit[3];
    vhat_c(pointOfInterestProj, pointOfInterestProjUnit);

    // Scale the vector to within a pixel of the azimuth's origin point.
    // Get pixel scale in km/pixel and divide by 2 to insure that we stay within
    // a pixel of the origin point
    double scale = (PixelResolution() / 1000.0) / 2.0;
    SpiceDouble pointOfInterestProjUnitScaled[3];
    vscl_c(scale, pointOfInterestProjUnit, pointOfInterestProjUnitScaled);

    // Compute the adjusted point of interest vector from the body center. This point
    // will be within a pixel of the origin and in the same direction as the requested
    // raw point of interest vector
    SpiceDouble adjustedPointOfInterestFromBodyCenter[3];
    vadd_c(azimuthOrigin, pointOfInterestProjUnitScaled, adjustedPointOfInterestFromBodyCenter);

    // Get the origin image coordinate
    double azimuthOriginSample = Sample();
    double azimuthOriginLine = Line();

    // Convert the point to a lat/lon and find out its image coordinate
    double adjustedPointOfInterestRad, adjustedPointOfInterestLon, adjustedPointOfInterestLat;
    reclat_c(adjustedPointOfInterestFromBodyCenter, 
             &adjustedPointOfInterestRad, 
             &adjustedPointOfInterestLon, 
             &adjustedPointOfInterestLat);
    adjustedPointOfInterestLat = adjustedPointOfInterestLat * 180.0 / PI;
    adjustedPointOfInterestLon = adjustedPointOfInterestLon * 180.0 / PI;
    if (adjustedPointOfInterestLon < 0) adjustedPointOfInterestLon += 360.0;

    // Use the radius of the azimuth's origin point
    // (rather than that of the adjusted point of interest location)
    // to avoid the effects of topography on the calculation
    bool success = SetUniversalGround(adjustedPointOfInterestLat, 
                                      adjustedPointOfInterestLon, 
                                      originRadius.meters());
    if (!success) {
      // if the adjusted lat/lon values for the point of origin fail to be set, we can not compute
      // an azimuth. reset to the original sample/line and return null.
      SetImage(azimuthOriginSample, azimuthOriginLine);
      return Isis::Null;
    }

    double adjustedPointOfInterestSample = Sample();
    double adjustedPointOfInterestLine   = Line();

    // TODO:  Write PushState and PopState method to ensure the
    // internals of the class are set based on SetImage or SetGround

    // We now have the information needed to calculate an arctangent
    // 
    // 
    //         point of interest
    //                  |\       |
    //                  | \      |
    //                  |  \     |                      tan(A) = (delta line) / (delta sample)
    //    delta line    |   \    |                      A  = arctan( (delta line) / (delta sample) )
    //                  |    \   |
    //                  |     \  |
    //                  |      \ |
    //       ___________|_____A_\|_______________
    //              delta sample |origin point
    //                           |
    //                           |
    //                           |
    //                           |
    //                           |
    //
    // in this example, the azimuth is the angle indicated by the A plus 180 degrees, since we begin
    // the angle rotation at the positive x-axis
    // This quadrant issue (i.e.need to add 180 degrees) is handled by the atan2 program
    // 
    double deltaSample = adjustedPointOfInterestSample - azimuthOriginSample;
    double deltaLine = adjustedPointOfInterestLine - azimuthOriginLine;
    
    // Compute the angle; the azimuth is the arctangent of the line difference divided by the
    // sample difference; the atan2 function is used because it determines which quadrant we
    // are in based on the sign of the 2 arguments; the arctangent is measured in a positive
    // clockwise direction because the lines in the image increase downward; the arctangent
    // uses the 3 o'clock axis (positive sample direction) as its reference line (line of
    // zero degrees); a good place to read about the atan2 function is at
    // http://en.wikipedia.org/wiki/Atan2
    double azimuth = 0.0;
    if (deltaSample != 0.0 || deltaLine != 0.0) {
      azimuth = atan2(deltaLine, deltaSample);
      azimuth *= 180.0 / PI;
    }

    // Azimuth is limited to the range of 0 to 360
    if (azimuth < 0.0) azimuth += 360.0;
    if (azimuth > 360.0) azimuth -= 360.0;

    NaifStatus::CheckErrors();

    // computed is true if the sample/line or lat/lon were reset in this method
    // to find the location of the point of interest
    // If so, reset "state" of camera to the original sample/line
    if (computed) {
      SetImage(azimuthOriginSample, azimuthOriginLine);
    }
    else {
      p_pointComputed = false;
    }

    return azimuth;
  }

  double Camera::OffNadirAngle() {
    NaifStatus::CheckErrors();

    // Get the xyz coordinates for the spacecraft and point we are interested in
    double coord[3], spCoord[3];
    Coordinate(coord);
    instrumentPosition(spCoord);

    // Get the angle between the 2 points and convert to degrees
    double a = vsep_c(coord, spCoord) * 180.0 / PI;
    double b = 180.0 - EmissionAngle();

    // The three angles in a triangle must add up to 180 degrees
    double c = 180.0 - (a + b);

    NaifStatus::CheckErrors();

    return c;
  }

  double Camera::GroundAzimuth(double glat, double glon,
                               double slat, double slon) {
    double a;
    double b;
    if (glat >= 0.0) {
      a = (90.0 - slat) * PI / 180.0;
      b = (90.0 - glat) * PI / 180.0;
    } 
    else {
      a = (90.0 + slat) * PI / 180.0;
      b = (90.0 + glat) * PI / 180.0;
    }

    double cslon = slon;
    double cglon = glon;
    if (cslon > cglon) {
      if ((cslon-cglon) > 180.0) {
        while ((cslon-cglon) > 180.0) cslon = cslon - 360.0;
      }
    }
    if (cglon > cslon) {
      if ((cglon-cslon) > 180.0) {
        while ((cglon-cslon) > 180.0) cglon = cglon - 360.0;
      }
    }
 
    // Which quadrant are we in?
    int quad;
    if (slat > glat) {
      if (cslon > cglon) {
        quad = 1;
      } 
      else if (cslon < cglon) {
        quad = 2;
      }
      else {
        quad = 1;
      }
    }
    else if (slat < glat) {
      if (cslon > cglon) {
        quad = 4;
      } 
      else if (cslon < cglon) {
        quad = 3;
      } 
      else {
        quad = 4;
      }
    } 
    else {
      if (cslon > cglon) {
        quad = 1;
      }
      else if (cslon < cglon) {
        quad = 2;
      }
      else {
        return 0.0;
      }
    }

    double C = (cglon - cslon) * PI / 180.0;
    if (C < 0) C = -C;
    double c = acos(cos(a)*cos(b) + sin(a)*sin(b)*cos(C));
    double azimuth = 0.0;
    if (sin(b) == 0.0 || sin(c) == 0.0) {
      return azimuth;
    }
    double A = acos((cos(a) - cos(b)*cos(c))/(sin(b)*sin(c))) * 180.0 / PI;
    //double B = acos((cos(b) - cos(c)*cos(a))/(sin(c)*sin(a))) * 180.0 / PI;
    if (glat >= 0.0) {
      if (quad == 1 || quad == 4) {
        azimuth = A;
      }
      else if (quad == 2 || quad == 3) {
        azimuth = 360.0 - A;
      }
    }
    else {
      if (quad == 1 || quad == 4) {
        azimuth = 180.0 - A;
      }
      else if (quad == 2 || quad == 3) {
        azimuth = 180.0 + A;
      }
    }

    return azimuth;
  }


  void Camera::SetDistortionMap(CameraDistortionMap *map) {
    if (p_distortionMap) {
      delete p_distortionMap;
    }

    p_distortionMap = map;
  }

  void Camera::SetFocalPlaneMap(CameraFocalPlaneMap *map) {
    if (p_focalPlaneMap) {
      delete p_focalPlaneMap;
    }

    p_focalPlaneMap = map;
  }

  void Camera::SetDetectorMap(CameraDetectorMap *map) {
    if (p_detectorMap) {
      delete p_detectorMap;
    }

    p_detectorMap = map;
  }

  void Camera::SetGroundMap(CameraGroundMap *map) {
    if (p_groundMap) {
      delete p_groundMap;
    }

    p_groundMap = map;
  }

  void Camera::SetSkyMap(CameraSkyMap *map) {
    if (p_skyMap) {
      delete p_skyMap;
    }

    p_skyMap = map;
  }

  void Camera::LoadCache() {
    // We want to stay in unprojected space for this process
    bool projIgnored = p_ignoreProjection;
    p_ignoreProjection = true;

    // get the cache variables
    pair<double,double> ephemerisTimes = StartEndEphemerisTimes();
    int cacheSize = CacheSize(ephemerisTimes.first, ephemerisTimes.second);

    // Set a position in the image so that the PixelResolution can be calculated
    SetImage(p_alphaCube->BetaSamples() / 2, p_alphaCube->BetaLines() / 2);
    double tol = PixelResolution() / 100.; //meters/pix/100.

    if (tol < 0.0) {
      // Alternative calculation of ground resolution of a pixel/100
      double altitudeMeters;
      if (target()->isSky()) {   // Use the unit sphere as the target
        altitudeMeters = 1.0;
      }
      else {
        altitudeMeters = SpacecraftAltitude() * 1000.;
      }
      tol = PixelPitch() * altitudeMeters / FocalLength() / 100.;
    }

    p_ignoreProjection = projIgnored;

    Spice::createCache(ephemerisTimes.first, ephemerisTimes.second,
                       cacheSize, tol);

    setTime(ephemerisTimes.first);

    // Reset to band 1
    SetBand(1);

    return;
  }


  pair<double, double> Camera::StartEndEphemerisTimes() {
    pair<double,double> ephemerisTimes;
    double startTime = -DBL_MAX;
    double endTime = -DBL_MAX;

    for (int band = 1; band <= Bands(); band++) {
      SetBand(band);
      SetImage(0.5, 0.5);
      double etStart = time().Et();
      SetImage(p_alphaCube->BetaSamples() + 0.5,
               p_alphaCube->BetaLines() + 0.5);
      double etEnd = time().Et();
      if (band == 1) {
        startTime = min(etStart, etEnd);
        endTime = max(etStart, etEnd);
      }
      startTime = min(startTime, min(etStart, etEnd));
      endTime = max(endTime, max(etStart, etEnd));
    }
    if (startTime == -DBL_MAX || endTime == -DBL_MAX) {
      string msg = "Unable to find time range for the spice kernels";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    ephemerisTimes.first = startTime;
    ephemerisTimes.second = endTime;
    return ephemerisTimes;
  }

  int Camera::CacheSize(double startTime, double endTime) {
    int cacheSize;
    // BetaLines() + 1 so we get at least 2 points for interpolation
    cacheSize = p_alphaCube->BetaLines() + 1;
    if (startTime == endTime) {
      cacheSize = 1;
    }
    return cacheSize;
  }


  void Camera::SetGeometricTilingHint(int startSize, int endSize) {
    // verify the start size is a multiple of 2 greater than 2
    int powerOf2 = 2;

    // No hint if 2's are passed in
    if (startSize == 2 && endSize == 2) {
      p_geometricTilingStartSize = 2;
      p_geometricTilingEndSize = 2;
      return;
    }

    if (endSize > startSize) {
      IString message = "Camera::SetGeometricTilingHint End size must be smaller than the start size";
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    if (startSize < 4) {
      IString message = "Camera::SetGeometricTilingHint Start size must be at least 4";
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    bool foundEnd = false;
    while (powerOf2 > 0 && startSize != powerOf2) {
      if (powerOf2 == endSize) foundEnd = true;
      powerOf2 *= 2;
    }

    // Didnt find a solution, the integer became negative first, must not be
    //   a power of 2
    if (powerOf2 < 0) {
      IString message = "Camera::SetGeometricTilingHint Start size must be a power of 2";
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    if (!foundEnd) {
      IString message = "Camera::SetGeometricTilingHint End size must be a power of 2 less than the start size, but greater than 2";
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    p_geometricTilingStartSize = startSize;
    p_geometricTilingEndSize = endSize;
  }

  void Camera::GetGeometricTilingHint(int &startSize, int &endSize) {
    startSize = p_geometricTilingStartSize;
    endSize = p_geometricTilingEndSize;
  }


  bool Camera::InCube() {
    if (Sample() < 0.5 || Line() < 0.5) {
      return false;
    }

    if (Sample() > Samples() + 0.5 || Line() > Lines() + 0.5) {
      return false;
    }

    return true;
  }

  bool Camera::HasProjection() {
    return p_projection != 0;
  }

  bool Camera::IsBandIndependent() {
    return true;
  }

  int Camera::ReferenceBand() const {
    return p_referenceBand;
  }

  bool Camera::HasReferenceBand() const {
    return p_referenceBand != 0;
  }

  void Camera::SetBand(const int band) {
    p_childBand = band;
  }

   double Camera::Sample() {
    return p_childSample;
  }

   int Camera::Band() {
    return p_childBand;
  }

   double Camera::Line() {
    return p_childLine;
  }
  double Camera::resolution() {
    return PixelResolution();
  }
   double Camera::FocalLength() const {
    return p_focalLength;
  }

   double Camera::PixelPitch() const {
    return p_pixelPitch;
  }

   QList<QPointF> Camera::PixelIfovOffsets() {

     QList<QPointF> offsets;
     offsets.append(QPointF(-PixelPitch() * DetectorMap()->SampleScaleFactor() / 2.0,
                            -PixelPitch() * DetectorMap()->LineScaleFactor() / 2.0));
     offsets.append(QPointF(PixelPitch() * DetectorMap()->SampleScaleFactor() / 2.0,
                            -PixelPitch() * DetectorMap()->LineScaleFactor() / 2.0));
     offsets.append(QPointF(PixelPitch() * DetectorMap()->SampleScaleFactor() / 2.0,
                            PixelPitch() * DetectorMap()->LineScaleFactor() / 2.0));
     offsets.append(QPointF(-PixelPitch() * DetectorMap()->SampleScaleFactor() / 2.0,
                            PixelPitch() * DetectorMap()->LineScaleFactor() / 2.0));

     return offsets;
   }


   int Camera::Samples() const {
    return p_samples;
  }

   int Camera::Lines() const {
    return p_lines;
  }

   int Camera::Bands() const {
    return p_bands;
  }

   int Camera::ParentLines() const {
    return p_alphaCube->AlphaLines();
  }

   int Camera::ParentSamples() const {
    return p_alphaCube->AlphaSamples();
  }
  CameraDistortionMap *Camera::DistortionMap() {
    return p_distortionMap;
  }

  CameraFocalPlaneMap *Camera::FocalPlaneMap() {
    return p_focalPlaneMap;
  }

  CameraDetectorMap *Camera::DetectorMap() {
    return p_detectorMap;
  }

  CameraGroundMap *Camera::GroundMap() {
    return p_groundMap;
  }

  CameraSkyMap *Camera::SkyMap() {
    return p_skyMap;
  }

  void Camera::IgnoreProjection(bool ignore) {
    p_ignoreProjection = ignore;
  }
  int Camera::SpkTargetId() const { 
    return (naifSpkCode());
  }

  int Camera::SpkCenterId() const { 
    return (naifBodyCode());
  }

  void Camera::SetFocalLength(double v) {
    p_focalLength = v;
  }

  void Camera::SetPixelPitch(double v) {
    p_pixelPitch = v;
  }
// end namespace isis
}