text: #include "Isis.h"
  
  #include <cmath>
  
  #include "Camera.h"
  #include "Cube.h"
  #include "Distance.h"
  #include "FileList.h"
  #include "Process.h"
  #include "Pvl.h"
  #include "Statistics.h"
  #include "Target.h"
  
  using namespace std;
  using namespace Isis;
  
  template <typename T> inline T MIN(const T &A, const T &B) {
    if(A < B) {
      return (A);
    }
    else         {
      return (B);
    }
  }
  
  template <typename T> inline T MAX(const T &A, const T &B) {
    if(A > B) {
      return (A);
    }
    else         {
      return (B);
    }
  }
  
  inline double SetFloor(double value, const int precision) {
    double scale = pow(10.0, precision);
    value = floor(value * scale) / scale;
    return (value);
  }
  
  inline double SetRound(double value, const int precision) {
    double scale = pow(10.0, precision);
    value = rint(value * scale) / scale;
    return (value);
  }
  
  inline double SetCeil(double value, const int precision) {
    double scale = pow(10.0, precision);
    value = ceil(value * scale) / scale;
    return (value);
  }
  
  inline double Scale(const double pixres, const double polarRadius,
                      const double equiRadius, const double trueLat = 0.0) {
    double lat = trueLat * Isis::PI / 180.0;
    double a = polarRadius * cos(lat);
    double b = equiRadius * sin(lat);
    double localRadius = equiRadius * polarRadius / sqrt(a * a + b * b);
    return (localRadius / pixres * pi_c() / 180.0);
  }
  
  void IsisMain() {
    Process p;
  
    // Get the list of names of input CCD cubes to stitch together
    FileList flist;
    UserInterface &ui = Application::GetUserInterface();
    flist.read(ui.GetFileName("FROMLIST"));
    if(flist.size() < 1) {
      QString msg = "The list file[" + ui.GetFileName("FROMLIST") +
                    " does not contain any filenames";
      throw IException(IException::User, msg, _FILEINFO_);
    }
  
    QString projection("Equirectangular");
    if(ui.WasEntered("MAP")) {
      Pvl mapfile(ui.GetFileName("MAP"));
      projection = (QString) mapfile.findGroup("Mapping")["ProjectionName"];
    }
  
    if(ui.WasEntered("PROJECTION")) {
      projection = ui.GetString("PROJECTION");
    }
  
    // Gather other user inputs to projection
    QString lattype = ui.GetString("LATTYPE");
    QString londir  = ui.GetString("LONDIR");
    QString londom  = ui.GetString("LONDOM");
    int digits = ui.GetInteger("PRECISION");
  
    // Fix them for mapping group
    lattype = (lattype == "PLANETOCENTRIC") ? "Planetocentric" : "Planetographic";
    londir = (londir == "POSITIVEEAST") ? "PositiveEast" : "PositiveWest";
  
    Progress prog;
    prog.SetMaximumSteps(flist.size());
    prog.CheckStatus();
  
    Statistics scaleStat;
  
    
    Statistics obliqueScaleStat;
  
    Statistics longitudeStat;
    Statistics latitudeStat;
    Statistics equiRadStat;
    Statistics poleRadStat;
    PvlObject fileset("FileSet");
  
    // Save major equitorial and polar radii for last occuring
    double eqRad;
    double poleRad;
  
    QString target("Unknown");
    for(int i = 0 ; i < flist.size() ; i++) {
      try {
        // Set the input image, get the camera model, and a basic mapping
        // group
        Cube cube;
        cube.open(flist[i].toString());
  
        int lines = cube.lineCount();
        int samples = cube.sampleCount();
  
  
        PvlObject fmap("File");
        fmap += PvlKeyword("Name", flist[i].toString());
        fmap += PvlKeyword("Lines", toString(lines));
        fmap += PvlKeyword("Samples", toString(samples));
  
        Camera *cam = cube.camera();
        Pvl mapping;
        cam->BasicMapping(mapping);
        PvlGroup &mapgrp = mapping.findGroup("Mapping");
        mapgrp.addKeyword(PvlKeyword("ProjectionName", projection), Pvl::Replace);
        mapgrp.addKeyword(PvlKeyword("LatitudeType", lattype), Pvl::Replace);
        mapgrp.addKeyword(PvlKeyword("LongitudeDirection", londir), Pvl::Replace);
        mapgrp.addKeyword(PvlKeyword("LongitudeDomain", londom), Pvl::Replace);
  
        // Get the radii
        Distance radii[3];
        cam->radii(radii);
  
        eqRad   = radii[0].meters();
        poleRad = radii[2].meters();
  
        target = cam->target()->name();
        equiRadStat.AddData(&eqRad, 1);
        poleRadStat.AddData(&poleRad, 1);
  
        // Get resolution
        double lowres = cam->LowestImageResolution();
        double hires = cam->HighestImageResolution();
  
        
        double lowObliqueRes = cam->LowestObliqueImageResolution();
  
  
        
        double hiObliqueRes= cam->HighestObliqueImageResolution();
  
  
        scaleStat.AddData(&hires, 1);
        scaleStat.AddData(&lowres, 1);
  
        
        obliqueScaleStat.AddData(&hiObliqueRes,1);
        obliqueScaleStat.AddData(&lowObliqueRes,1);
  
  
        double pixres = (lowres + hires) / 2.0;
  
        
        //double obliquePixRes = (lowObliqueRes+hiObliqueRes)/2.0;
  
        double scale = Scale(pixres, poleRad, eqRad);
  
        
        //double obliqueScale = Scale(obliquePixRes,poleRad,eqRad);
  
  
        mapgrp.addKeyword(PvlKeyword("PixelResolution", toString(pixres)), Pvl::Replace);
  
        
        //mapgrp.addKeyword(PvlKeyword("ObliquePixelResolution", toString(obliquePixRes)),
        //                  Pvl::Replace);
  
        mapgrp.addKeyword(PvlKeyword("Scale", toString(scale), "pixels/degree"), Pvl::Replace);
  
        
        //mapgrp.addKeyword(PvlKeyword("ObliqueScale", toString(obliqueScale), "pixels/degree"),
        //                  Pvl::Replace);
        mapgrp += PvlKeyword("MinPixelResolution", toString(lowres), "meters/pixel");
        mapgrp += PvlKeyword("MaxPixelResolution", toString(hires), "meters/pixel");
  
        
        mapgrp += PvlKeyword("MinObliquePixelResolution", toString(lowObliqueRes), "meters/pixel");
        mapgrp += PvlKeyword("MaxObliquePixelResolution", toString(hiObliqueRes), "meters/pixel");
  
        // Get the universal ground range
        double minlat, maxlat, minlon, maxlon;
        cam->GroundRange(minlat, maxlat, minlon, maxlon, mapping);
        mapgrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace);
        mapgrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace);
        mapgrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace);
        mapgrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace);
  
        fmap.addGroup(mapgrp);
        fileset.addObject(fmap);
  
        longitudeStat.AddData(&minlon, 1);
        longitudeStat.AddData(&maxlon, 1);
        latitudeStat.AddData(&minlat, 1);
        latitudeStat.AddData(&maxlat, 1);
      }
      catch(IException &ie) {
        QString mess = "Problems with file " + flist[i].toString() + "\n" +
                       ie.what();
        throw IException(IException::User, mess, _FILEINFO_);
      }
  
      p.ClearInputCubes();
      prog.CheckStatus();
    }
  
  //  Construct the output mapping group with statistics
    PvlGroup mapping("Mapping");
    double avgPixRes( (scaleStat.Minimum() + scaleStat.Maximum() ) / 2.0);
  
    
    //double avgObliquePixRes( (obliqueScaleStat.Minimum() + obliqueScaleStat.Maximum() ) / 2.0);
  
    double avgLat((latitudeStat.Minimum() + latitudeStat.Maximum()) / 2.0);
    double avgLon((longitudeStat.Minimum() + longitudeStat.Maximum()) / 2.0);
    double avgEqRad((equiRadStat.Minimum() + equiRadStat.Maximum()) / 2.0);
    double avgPoleRad((poleRadStat.Minimum() + poleRadStat.Maximum()) / 2.0);
    double scale  = Scale(avgPixRes, avgPoleRad, avgEqRad);
    //double obliqueScale = Scale(avgObliquePixRes,avgPoleRad,avgEqRad);
  
    mapping += PvlKeyword("ProjectionName", projection);
    mapping += PvlKeyword("TargetName", target);
    mapping += PvlKeyword("EquatorialRadius", toString(eqRad), "meters");
    mapping += PvlKeyword("PolarRadius", toString(poleRad), "meters");
    mapping += PvlKeyword("LatitudeType", lattype);
    mapping += PvlKeyword("LongitudeDirection", londir);
    mapping += PvlKeyword("LongitudeDomain", londom);
    mapping += PvlKeyword("PixelResolution", toString(SetRound(avgPixRes, digits)), "meters/pixel");
  
    
    //mapping += PvlKeyword("ObliquePixelResolution", toString(SetRound(avgObliquePixRes, digits)),
    //                      "meters/pixel");
  
  
    mapping += PvlKeyword("Scale", toString(SetRound(scale, digits)), "pixels/degree");
  
    
    //mapping += PvlKeyword("ObliqueScale", toString(SetRound(obliqueScale, digits)), "pixels/degree");
  
  
    mapping += PvlKeyword("MinPixelResolution", toString(scaleStat.Minimum()), "meters/pixel");
    mapping += PvlKeyword("MaxPixelResolution", toString(scaleStat.Maximum()), "meters/pixel");
  
    
    mapping += PvlKeyword("MinObliquePixelResolution", toString(obliqueScaleStat.Minimum()),
                          "meters/pixel");
    mapping += PvlKeyword("MaxObliquePixelResolution", toString(obliqueScaleStat.Maximum()),
                          "meters/pixel");
  
    mapping += PvlKeyword("CenterLongitude", toString(SetRound(avgLon, digits)));
    mapping += PvlKeyword("CenterLatitude",  toString(SetRound(avgLat, digits)));
    mapping += PvlKeyword("MinimumLatitude", toString(MAX(SetFloor(latitudeStat.Minimum(),
                                                                   digits), -90.0)));
    mapping += PvlKeyword("MaximumLatitude", toString(MIN(SetCeil(latitudeStat.Maximum(),
                                                                   digits), 90.0)));
    mapping += PvlKeyword("MinimumLongitude", toString(MAX(SetFloor(longitudeStat.Minimum(),
                                                                    digits), -180.0)));
    mapping += PvlKeyword("MaximumLongitude", toString(MIN(SetCeil(longitudeStat.Maximum(),
                                                                   digits), 360.0)));
  
    PvlKeyword clat("PreciseCenterLongitude", toString(avgLon));
    clat.addComment("Actual Parameters without precision applied");
    mapping += clat;
    mapping += PvlKeyword("PreciseCenterLatitude",  toString(avgLat));
    mapping += PvlKeyword("PreciseMinimumLatitude", toString(latitudeStat.Minimum()));
    mapping += PvlKeyword("PreciseMaximumLatitude", toString(latitudeStat.Maximum()));
    mapping += PvlKeyword("PreciseMinimumLongitude", toString(longitudeStat.Minimum()));
    mapping += PvlKeyword("PreciseMaximumLongitude", toString(longitudeStat.Maximum()));
  
  
    Application::Log(mapping);
  
    // Write the output file if requested
    if(ui.WasEntered("TO")) {
      Pvl temp;
      temp.addGroup(mapping);
      temp.write(ui.GetFileName("TO", "map"));
    }
  
    if(ui.WasEntered("LOG")) {
      Pvl temp;
      temp.addObject(fileset);
      temp.write(ui.GetFileName("LOG", "log"));
    }
  
    p.EndProcess();
  }
  