Unverified Commit e64f1566 authored by Alessandro Frigeri's avatar Alessandro Frigeri Committed by GitHub
Browse files

Refactor of the pixel2map application and add geospatial vector file output (#5606)

* added README.md

* added pixel2map '23 mods

* updated the 2023 code to the current dev branch

* Started the refactoring

* Added app functional test

* updated files for refactoring

* removed call to UserInterface in new cpp code

* refactored started working

* First functional test working

* Removed unused headers

* code reviewed with Kelvin

* fixed Functional Test - started UnitTest

* added the unit test cpp file in the tests folder

* added new test, fixed code in case of 360 crossing

* cleanup for PR

* updated CHANGELOG and docs

* updated institution info

* removed commented-out lines in header

* removed README.md

* pixel2map definition re-inserted into Isis namespace (which was deleted during code cleanup)

* change on testCube path reverted
parent 0dff5234
Loading
Loading
Loading
Loading
+5 −0
Original line number Original line Diff line number Diff line
@@ -166,6 +166,11 @@
      "affiliation": "United States Geological Survey, Astro Geology Science Center",
      "affiliation": "United States Geological Survey, Astro Geology Science Center",
      "name": "Fergason, Robin"
      "name": "Fergason, Robin"
    },
    },
    {
      "affiliation": "Italian National Institute for Astrophysics (INAF), Istituto di Astrofisica e Planetologia Spaziali (IAPS), Rome, Italy",
      "name": "Frigeri, Alessandro",
      "orcid": "0000-0002-9140-3977"   
    },
    {
    {
      "name": "Gaddis, Lisa"
      "name": "Gaddis, Lisa"
    },
    },
+9 −1
Original line number Original line Diff line number Diff line
@@ -35,8 +35,16 @@ release.


## [Unreleased]
## [Unreleased]


### Fixed
### Added
- Added TOVECT output parameter which generate a geospatial CSV file with a VRT metadata sidecar file [#5571](https://github.com/DOI-USGS/ISIS3/issues/5571)  
- Added Vectorize to ProcessGroundPolygon library
- Added gtest files for the app and unit test 


### Changed
- Refactored the pixel2map app
- Updated pixel2map documentation

### Fixed
- Fixed a bug in kaguyasp2isis that doesn't work for data with a detached label.
- Fixed a bug in kaguyasp2isis that doesn't work for data with a detached label.


## [8.3.0] - 2024-08-16
## [8.3.0] - 2024-08-16
+16 −385
Original line number Original line Diff line number Diff line
#define GUIHELPERS
#define GUIHELPERS

#include "Isis.h"
#include "Isis.h"


#include <QDebug>
#include "Application.h"
#include <QList>
#include <QPointF>
#include <QString>

#include "Camera.h"
#include "Cube.h"
#include "CubeAttribute.h"
#include "FileList.h"
#include "IException.h"
#include "PixelFOV.h"
#include "ProcessByBrick.h"
#include "ProcessGroundPolygons.h"
#include "ProcessRubberSheet.h"
#include "ProjectionFactory.h"
#include "Pvl.h"
#include "PvlGroup.h"
#include "Target.h"


#include "pixel2map.h"
#include "pixel2map.h"


#include "UserInterface.h"

using namespace std;
using namespace std;
using namespace Isis;
using namespace Isis;


void PrintMap();
void PrintMap();
void rasterizePixel(Isis::Buffer &in);


map <QString, void *> GuiHelpers() {
// Check for docs for GUI helpers -> KEEP and move it to pixel2map.cpp 
  std::map <QString, void *> GuiHelpers() {
  map <QString, void *> helper;
  map <QString, void *> helper;
  helper ["PrintMap"] = (void *) PrintMap;
  helper ["PrintMap"] = (void *) PrintMap;
  return helper;
  return helper;
}
}


// Global variables
ProcessGroundPolygons g_processGroundPolygons;
Camera *g_incam;
int g_numIFOVs = 0;


void IsisMain() {
void IsisMain() {

  g_incam = NULL;

  // Get the map projection file provided by the user
  UserInterface &ui = Application::GetUserInterface();
  UserInterface &ui = Application::GetUserInterface();
  Pvl userMap;
  Pvl appLog;
  userMap.read(ui.GetFileName("MAP"));
  pixel2map(ui, &appLog);
  PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse);

  FileList list;
  if (ui.GetString("FROMTYPE") == "FROM") {
    // GetAsString will capture the entire string, including attributes
    list.push_back(FileName(ui.GetAsString("FROM")));
  }
  else {
    try {
      list.read(ui.GetFileName("FROMLIST"));
    }
    catch (IException &e) {
      throw IException(e);
    }
  }

  if (ui.GetString("FOVRANGE") == "INSTANTANEOUS") {
    g_numIFOVs = 1;
  }
  else {
    g_numIFOVs = ui.GetInteger("NUMIFOV");
  }

  double newminlat, newmaxlat, newminlon, newmaxlon;
  double minlat = 90;
  double maxlat = -90;
  double minlon = 360.0;
  double maxlon = 0;
  PvlGroup camGrp;
  QString lastBandString;

  // Get the combined lat/lon range for all input cubes
  int bands = 1;
  for (int i = 0; i < list.size(); i++) {

    // Open the input cube and get the camera
    CubeAttributeInput atts0(list[i]);
    Cube icube;
    if(atts0.bands().size() != 0) {
      vector<QString> lame = atts0.bands();
      icube.setVirtualBands(lame);
    }
    icube.open( list[i].toString() );
    bands = icube.bandCount();
    g_incam = icube.camera();

    // Make sure it is not the sky
    if (g_incam->target()->isSky()) {
      QString msg = "The image [" + list[i].toString() +
                    "] is targeting the sky, use skymap instead.";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    // Make sure all the bands for all the files match
    if (i >= 1 && atts0.bandsString() != lastBandString) {
      QString msg = "The Band numbers for all the files do not match.";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    else {
      lastBandString = atts0.bandsString();
    }

    // Get the mapping group
    Pvl camMap;
    g_incam->BasicMapping(camMap);
    camGrp = camMap.findGroup("Mapping");

    g_incam->GroundRange(newminlat, newmaxlat, newminlon, newmaxlon, userMap);
    //set min lat/lon
    if (newminlat < minlat) {
      minlat = newminlat;
    }
    if (newminlon < minlon) {
      minlon = newminlon;
    }

    //set max lat/lon
    if (newmaxlat > maxlat) {
      maxlat = newmaxlat;
    }
    if (newmaxlon > maxlon) {
      maxlon = newmaxlon;
    }

    // The camera will be deleted when the Cube is deleted so NULL g_incam
    g_incam = NULL;
  } //end for list.size

  camGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace);
  camGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace);
  camGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace);
  camGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace);


  // We want to delete the keywords we just added if the user wants the range
  // out of the mapfile, otherwise they will replace any keywords not in the
  // mapfile
  if (ui.GetString("DEFAULTRANGE") == "MAP") {
    camGrp.deleteKeyword("MinimumLatitude");
    camGrp.deleteKeyword("MaximumLatitude");
    camGrp.deleteKeyword("MinimumLongitude");
    camGrp.deleteKeyword("MaximumLongitude");
  }

  // Otherwise, remove the keywords from the map file so the camera keywords
  // will be propogated correctly
  else {
    while (userGrp.hasKeyword("MinimumLatitude")) {
      userGrp.deleteKeyword("MinimumLatitude");
    }
    while (userGrp.hasKeyword("MinimumLongitude")) {
      userGrp.deleteKeyword("MinimumLongitude");
    }
    while (userGrp.hasKeyword("MaximumLatitude")) {
      userGrp.deleteKeyword("MaximumLatitude");
    }
    while (userGrp.hasKeyword("MaximumLongitude")) {
      userGrp.deleteKeyword("MaximumLongitude");
    }
  }

  // If the user decided to enter a ground range then override
  if (ui.WasEntered("MINLON")) {
    userGrp.addKeyword(PvlKeyword("MinimumLongitude",
                                  toString(ui.GetDouble("MINLON"))), Pvl::Replace);
  }

  if (ui.WasEntered("MAXLON")) {
    userGrp.addKeyword(PvlKeyword("MaximumLongitude",
                                  toString(ui.GetDouble("MAXLON"))), Pvl::Replace);
  }

  if (ui.WasEntered("MINLAT")) {
    userGrp.addKeyword(PvlKeyword("MinimumLatitude",
                                  toString(ui.GetDouble("MINLAT"))), Pvl::Replace);
  }

  if (ui.WasEntered("MAXLAT")) {
    userGrp.addKeyword(PvlKeyword("MaximumLatitude",
                                  toString(ui.GetDouble("MAXLAT"))), Pvl::Replace);
  }

  // If they want the res. from the mapfile, delete it from the camera so
  // nothing gets overriden
  if (ui.GetString("PIXRES") == "MAP") {
    camGrp.deleteKeyword("PixelResolution");
  }
  // Otherwise, delete any resolution keywords from the mapfile so the camera
  // info is propogated over
  else if (ui.GetString("PIXRES") == "CAMERA") {
    if (userGrp.hasKeyword("Scale")) {
      userGrp.deleteKeyword("Scale");
    }
    if (userGrp.hasKeyword("PixelResolution")) {
      userGrp.deleteKeyword("PixelResolution");
    }
  }
  
  
  // Copy any defaults that are not in the user map from the camera map file
  if( ui.WasEntered("TO") && ui.IsInteractive() ) {
  for (int k = 0; k < camGrp.keywords(); k++) {
    Application::GuiLog(appLog);
    if (!userGrp.hasKeyword(camGrp[k].name())) {
      userGrp += camGrp[k];
  }
  }
  }

  // If the user decided to enter a resolution then override
  if (ui.GetString("PIXRES") == "MPP") {
    userGrp.addKeyword(PvlKeyword("PixelResolution",
                                  toString(ui.GetDouble("RESOLUTION"))),
                       Pvl::Replace);
    if (userGrp.hasKeyword("Scale")) {
      userGrp.deleteKeyword("Scale");
    }
  }
  else if (ui.GetString("PIXRES") == "PPD") {
    userGrp.addKeyword(PvlKeyword("Scale",
                                  toString(ui.GetDouble("RESOLUTION"))),
                       Pvl::Replace);
    if (userGrp.hasKeyword("PixelResolution")) {
      userGrp.deleteKeyword("PixelResolution");
    }
  }

  // See if the user want us to handle the longitude seam
  if (ui.GetString("DEFAULTRANGE") == "CAMERA" || ui.GetString("DEFAULTRANGE") == "MINIMIZE") {

    // Open the last cube and use its camera
    CubeAttributeInput atts0( list.back() );
    Cube icube;
    if(atts0.bands().size() != 0) {
      vector<QString> lame = atts0.bands();
      icube.setVirtualBands(lame);
    }
    icube.open( list.back().toString() );
    g_incam = icube.camera();

    if (g_incam->IntersectsLongitudeDomain(userMap)) {
      if (ui.GetString("LONSEAM") == "AUTO") {
        if ((int) userGrp["LongitudeDomain"] == 360) {
          userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(180)),
                             Pvl::Replace);
          if (g_incam->IntersectsLongitudeDomain(userMap)) {
            // Its looks like a global image so switch back to the
            // users preference
            userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(360)),
                               Pvl::Replace);
          }
        }
        else {
          userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(360)),
                             Pvl::Replace);
          if (g_incam->IntersectsLongitudeDomain(userMap)) {
            // Its looks like a global image so switch back to the
            // users preference
            userGrp.addKeyword(PvlKeyword("LongitudeDomain", toString(180)),
                               Pvl::Replace);
          }
        }
        // Make the target info match the new longitude domain
        double minlat, maxlat, minlon, maxlon;
        g_incam->GroundRange(minlat, maxlat, minlon, maxlon, userMap);
        if (!ui.WasEntered("MINLAT")) {
          userGrp.addKeyword(PvlKeyword("MinimumLatitude", toString(minlat)), Pvl::Replace);
        }
        if (!ui.WasEntered("MAXLAT")) {
          userGrp.addKeyword(PvlKeyword("MaximumLatitude", toString(maxlat)), Pvl::Replace);
        }
        if (!ui.WasEntered("MINLON")) {
          userGrp.addKeyword(PvlKeyword("MinimumLongitude", toString(minlon)), Pvl::Replace);
        }
        if (!ui.WasEntered("MAXLON")) {
          userGrp.addKeyword(PvlKeyword("MaximumLongitude", toString(maxlon)), Pvl::Replace);
        }
      }

      else if (ui.GetString("LONSEAM") == "ERROR") {
        QString msg = "The image [" + ui.GetCubeName("FROM") + "] crosses the " +
                     "longitude seam";
        throw IException(IException::User, msg, _FILEINFO_);
      }
    }

    // The camera will be deleted when the Cube is deleted so NULL g_incam
    g_incam = NULL;
  }


  Pvl pvl;
  pvl.addGroup(userGrp);

  // If there is only one input cube, we need to attach an AlphaCube to the outputs.
  if (list.size() == 1) {
    Cube parent(list[0].toString());
    if (!parent.hasGroup("AlphaCube")) {
      PvlGroup alpha("AlphaCube");
      alpha += PvlKeyword("AlphaSamples", toString(parent.sampleCount()));
      alpha += PvlKeyword("AlphaLines", toString(parent.lineCount()));
      alpha += PvlKeyword("AlphaStartingSample", toString(0.5));
      alpha += PvlKeyword("AlphaStartingLine", toString(0.5));
      alpha += PvlKeyword("AlphaEndingSample", toString(parent.sampleCount() + 0.5));
      alpha += PvlKeyword("AlphaEndingLine", toString(parent.lineCount() + 0.5));
      alpha += PvlKeyword("BetaSamples", toString(parent.sampleCount()));
      alpha += PvlKeyword("BetaLines", toString(parent.lineCount()));
      pvl.addGroup(alpha);
    }
  }

  g_processGroundPolygons.SetStatCubes("TO", pvl, bands);

  bool useCenter = true;
  if (ui.GetString("METHOD") == "CENTER") {
    useCenter = true;
  }
  else if (ui.GetString("METHOD") == " WHOLEPIXEL") {
    useCenter = false;
  }
 
  g_processGroundPolygons.SetIntersectAlgorithm(useCenter);

  for (int f = 0; f < list.size(); f++) {

    Cube cube(list[f].toString(), "r");
    // Loop through the input cube and get the all pixels values for all bands
    ProcessByBrick processBrick;
    processBrick.Progress()->SetText("Working on file:  " + list[f].toString());
    processBrick.SetBrickSize(1, 1, bands);
    // Recall list[f] is a FileName, which stores the attributes
    CubeAttributeInput atts0(list[f]);
    Cube *icube = processBrick.SetInputCube(list[f].toString(), atts0, 0);
    g_incam = icube->camera();

    processBrick.StartProcess(rasterizePixel);
    processBrick.EndProcess();
  }
  
  // When there is only one input cube, we want to propagate IsisCube labels to output cubes
  if (list.size() == 1) {
    // Note that polygons and original labels are not propagated
    g_processGroundPolygons.PropagateLabels(list[0].toString());
    // Tell Process which tables we want to propagate
    QList<QString> tablesToPropagate;
    tablesToPropagate << "InstrumentPointing" << "InstrumentPosition" << "BodyRotation"
        << "SunPosition";
    g_processGroundPolygons.PropagateTables(list[0].toString(), tablesToPropagate);
  }
  g_processGroundPolygons.EndProcess();

  // WARNING: rasterizePixel() method alters the current state of the camera.
  // If any code is added after this point, you must call setImage to return
  // to original camera state before rasterization.
    
    
}
}


@@ -366,56 +36,17 @@ void IsisMain() {
  * Helper function to print out mapfile to session log
  * Helper function to print out mapfile to session log
  */
  */
void PrintMap() {
void PrintMap() {
//removed in the refactoring process 
UserInterface &ui = Application::GetUserInterface();
UserInterface &ui = Application::GetUserInterface();


  // Get mapping group from map file
  // Get mapping group from map file
  Pvl userMap;
  Pvl userMap(ui.GetFileName("MAP"));
  userMap.read(ui.GetFileName("MAP"));
  PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse);
  PvlGroup &userGrp = userMap.findGroup("Mapping", Pvl::Traverse);


  //Write map file out to the log
  //Write map file out to the log
//
  Isis::Application::GuiLog(userGrp);
  Isis::Application::GuiLog(userGrp);
}
} // PrintMap



/**
  * This method uses the ProcessGroundPolygons object to rasterize each 
  * pixel in the given buffer. 
  *  
  * @param in Input ProcessByBrick buffer. 
  */
void rasterizePixel(Isis::Buffer &in) {

  std::vector<double>lat, lon;
  std::vector<double>dns;

  for (int i = 0; i < in.size(); i++) {
    dns.push_back(in[i]);
  }
  
  
  int l = in.Line();
  int s = in.Sample();
  // TODO: This needs to be done for each band for band dependant instruments
  // Note: This can slow this down a lot

  // Get the IFOVs in lat/lon space
  PixelFOV fov;
  QList< QList< QPointF > > fovVertices = fov.latLonVertices(*g_incam, l, s, g_numIFOVs);

  // loop through each ifov list
  for (int ifov = 0; ifov < fovVertices.size(); ifov++) {
    // we need at least 3 vertices for a polygon
    if (fovVertices[ifov].size() > 3) {
      //  Get lat/lon for each vertex of the ifov
      for (int point = 0; point < fovVertices[ifov].size(); point++) {
        lat.push_back(fovVertices[ifov][point].x());
        lon.push_back(fovVertices[ifov][point].y());
      }
      // rasterize this ifov and clear vectors for next ifov
      g_processGroundPolygons.Rasterize(lat, lon, dns);
      lat.clear();
      lon.clear();
    }
  }
}
+500 −0

File added.

Preview size limit exceeded, changes collapsed.

+21 −18
Original line number Original line Diff line number Diff line
#ifndef pixel2map_h
#ifndef pixel2map_h
#define pixel2map_h
#define pixel2map_h


#include <QDebug>
#include <QList>
#include <QPointF>
#include <QString>

#include "Camera.h"
#include "Cube.h"
#include "FileList.h"
#include "PixelFOV.h"
#include "ProcessByBrick.h"
#include "ProcessGroundPolygons.h"
#include "Pvl.h"
#include "Target.h"

#include "Transform.h"
#include "Transform.h"


/**
#include "geos/geom/Geometry.h"
 * @author 2008-02-13 Stacy Alley
#include "geos/io/WKTWriter.h"
 *

 * @internal 
#include "UserInterface.h"
 *   @history 2013-07-30 Stuart C. Sides & Tracie Sucharski Renamed from vim2map
 */
class pixel2map : public Isis::Transform {
  private:
    Isis::Camera *p_incam;
    Isis::Projection *p_outmap;


  public:
    // constructor
    pixel2map(const int inputSamples, const int inputLines, Isis::Camera *incam,
             const int outputSamples, const int outputLines, Isis::Projection *outmap,
             bool trim);


    // destructor
namespace Isis {
    ~pixel2map() {};
	
	
extern void pixel2map(UserInterface &ui, Pvl *log=nullptr);


};
} 


#endif
#endif
Loading