Commit 806ba31d authored by Curtis Rose's avatar Curtis Rose
Browse files

Fixed an error where nocam2map would get into an infinite loop with cubes with...

Fixed an error where nocam2map would get into an infinite loop with cubes with a sample or line size less than 10. Fixes #2284.

git-svn-id: http://subversion.wr.usgs.gov/repos/prog/isis3/trunk@6787 41f8697f-d340-4b68-9986-7bafba869bb8
parent 7e9eae3a
Loading
Loading
Loading
Loading
+71 −15
Original line number Diff line number Diff line
@@ -21,9 +21,9 @@
#include "ProcessRubberSheet.h"
#include "ProjectionFactory.h"
#include "Statistics.h"
#include "Target.h"
#include "TextFile.h"
#include "TProjection.h"
#include "NaifStatus.h"

using namespace std;
using namespace Isis;
@@ -62,16 +62,15 @@ void IsisMain() {
  UserInterface &ui = Application::GetUserInterface();

  //Set the sample and line increments
  int sinc = (int)(inCube->sampleCount() * 0.10);
  float sinc = (inCube->sampleCount() * 0.10);
  if (ui.WasEntered("SINC")) {
    sinc = ui.GetInteger("SINC");
  }

  int linc = (int)(inCube->lineCount() * 0.10);
  float linc = (inCube->lineCount() * 0.10);
  if (ui.WasEntered("LINC")) {
    linc = ui.GetInteger("LINC");
  }

  //Set the degree of the polynomial to use in our functions
  int degree = ui.GetInteger("DEGREE");

@@ -83,6 +82,7 @@ void IsisMain() {
  LeastSquares sampSol(sampFunct);
  LeastSquares lineSol(lineFunct);


  //Setup the variables for solving the stereographic projection
  //x = cos(latitude) * sin(longitude - lon_center)
  //y = cos(lat_center) * sin(latitude) - sin(lat_center) * cos(latitude) * cos(longitude - lon_center)
@@ -90,15 +90,13 @@ void IsisMain() {
  //Get the center lat and long from the input cubes
  double lat_center = latCube->statistics()->Average() * PI / 180.0;
  double lon_center = lonCube->statistics()->Average() * PI / 180.0;


  /**
   * Loop through lines and samples projecting the latitude and longitude at those
   * points to stereographic x and y and adding these points to the LeastSquares
   * matrix.
   */
  for (int i = 1; i <= inCube->lineCount(); i += linc) {
    for (int j = 1; j <= inCube->sampleCount(); j += sinc) {
  for (float i = 1; i <= inCube->lineCount(); i += linc) {
    for (float j = 1; j <= inCube->sampleCount(); j += sinc) {
      latBrick.SetBasePosition(j, i, 1);
      latCube->read(latBrick);
      if(IsSpecial(latBrick.at(0))) continue;
@@ -134,8 +132,15 @@ void IsisMain() {
  }

  //Solve the least squares functions using QR Decomposition
  try {
    sampSol.Solve(LeastSquares::QRD);
    lineSol.Solve(LeastSquares::QRD);
  }
  catch (IException &e) {
    FileName inFile = ui.GetFileName("FROM");  
    QString msg = "Unable to calculate transformation of projection for [" + inFile.expanded() + "].";
    throw IException(e, IException::Unknown, msg, _FILEINFO_);
  }

  //If the user wants to save the residuals to a file, create a file and write
  //the column titles to it.
@@ -144,7 +149,6 @@ void IsisMain() {
    oFile.Open(ui.GetFileName("RESIDUALS"), "overwrite");
    oFile.PutLine("Sample,\tLine,\tX,\tY,\tSample Error,\tLine Error\n");
  }

  //Gather the statistics for the residuals from the least squares solutions
  Statistics sampErr;
  Statistics lineErr;
@@ -226,6 +230,8 @@ void IsisMain() {

    PvlKeyword equRadius;
    PvlKeyword polRadius;


    //If the user entered the equatorial and polar radii
    if (ui.WasEntered("EQURADIUS") && ui.WasEntered("POLRADIUS")) {
      equRadius = PvlKeyword("EquatorialRadius", toString(ui.GetDouble("EQURADIUS")));
@@ -233,10 +239,36 @@ void IsisMain() {
    }
    //Else read them from the pck
    else {
      PvlGroup radii = Target::radiiGroup(targetName[0]);
      equRadius = radii["EquatorialRadius"];
      polRadius = radii["PolarRadius"];
      FileName pckFile("$base/kernels/pck/pck?????.tpc");
      pckFile = pckFile.highestVersion();

      QString pckFileName = pckFile.expanded();

      NaifStatus::CheckErrors();
      furnsh_c(pckFileName.toAscii().data());

      QString target = targetName[0];
      SpiceInt code;
      SpiceBoolean found;

      bodn2c_c(target.toAscii().data(), &code, &found);
      NaifStatus::CheckErrors();

      if (!found) {
        QString msg = "Could not convert Target [" + target +
                     "] to NAIF code";
        throw IException(IException::Io, msg, _FILEINFO_);
      }

      SpiceInt n;
      SpiceDouble radii[3];

      bodvar_c(code, "RADII", &n, radii);

      equRadius = PvlKeyword("EquatorialRadius", toString(radii[0] * 1000));
      polRadius = PvlKeyword("PolarRadius", toString(radii[2] * 1000));
    }

    mapGrp.addKeyword(equRadius, Pvl::Replace);
    mapGrp.addKeyword(polRadius, Pvl::Replace);

@@ -877,6 +909,14 @@ void ComputeInputRange() {
    }
    //Else read them from the pck
    else {
      FileName pckFile("$base/kernels/pck/pck?????.tpc");
      pckFile = pckFile.highestVersion();

      QString pckFileName = pckFile.expanded();

      NaifStatus::CheckErrors();
      furnsh_c(pckFileName.toAscii().data());

      QString target;

      //If user entered target
@@ -890,9 +930,25 @@ void ComputeInputRange() {
        target = (QString)fromFile.findKeyword("TargetName", Pvl::Traverse);
      }

      PvlGroup radii = Target::radiiGroup(target);
      equRadius = double(radii["EquatorialRadius"]);
      polRadius = double(radii["PolarRadius"]);
      SpiceInt code;
      SpiceBoolean found;

      bodn2c_c(target.toAscii().data(), &code, &found);
      NaifStatus::CheckErrors();

      if (!found) {
        QString msg = "Could not convert Target [" + target +
                     "] to NAIF code";
        throw IException(IException::Io, msg, _FILEINFO_);
      }

      SpiceInt n;
      SpiceDouble radii[3];

      bodvar_c(code, "RADII", &n, radii);

      equRadius = radii[0] * 1000;
      polRadius = radii[2] * 1000;
    }

    if (isOcentric) {
+19 −12
Original line number Diff line number Diff line
@@ -27,13 +27,13 @@
    </p>

    <p>
      Therefore the root of the problem is finding the latitude/longitude in the 
      Therefore, the root of the problem is finding the latitude/longitude in the 
      two input cubes as it is a 2-dimensional search. We do this by first
      creating two polynomial functions of a user-specified order (DEGREE), one
      for sample and one for line. These functions, f(lat,lon) = sample and
      g(lat,lon) = line, are derived using a least-squares computation by
      collecting a sparse set of control points (line,samp)-to-(lat,lon). The
      sparseness is controlled by the SINC and LINC. Unfortunately the
      sparseness is controlled by the SINC and LINC. Unfortunately, the
      polynomial is not accurate enough to supply the mapping to sub-pixel
      accuracy. The level of errors between the control points and the derived
      equations can be written to the RESIDUAL file for examination.
@@ -97,6 +97,13 @@
      generated by the projection function.
    </p>

    <p>
      This application fails to produce correct output with smaller cubes sized 5x1, 1x5 or smaller. 
      The error message provided is: "**ERROR** Unable to solve-least squares using QR Decomposition.
      The upper triangular R matrix is not full rank in LeastSquares.cpp at 402."
      This issue occurs because there is not enough information provided by small cubes for the equations used.
    </p>

  </description>

  <category>
@@ -129,9 +136,9 @@
    <change name="Kristin Berry" date="2015-07-22">
      Added NaifStatus::CheckErrors() to see if any NAIF errors were signaled. References #2248
    </change>
    <change name="Jeannie Backer" date="2016-05-10">
      Replaced calls to NAIF routines (bodn2c and bodvar) with call to static method
      Isis::Target::radiiGroup. References #3934
    <change name="Curtis Rose" date="2016-05-16">
      Fixed an issue with small cubes (9x9 or less) causing application to become stuck in an infinite loop. References #2284.
      Also added to the description a message about why the application fails on small cubes (5x1, 1x5, or less).
    </change>
  </history>

@@ -500,7 +507,7 @@
          This parameter is used to specify how the default latitude/longitude
          ground range for the output map projected cube is obtained.  The
          ground range can be obtained from the camera or map file.  Note the
          user can overide the default using the MINLAT, MAXLAT, MINLON,
          user can override the default using the MINLAT, MAXLAT, MINLON,
          MAXLON parameters.  The purpose of the ground range is to define the
          coverage of the map projected cube.  Essentially, the ground range and
          pixel resolution are used to compute the size (samples and line) of the
@@ -529,7 +536,7 @@
           <option value="COMPUTE">
              <brief> Compute default range from input cube</brief>
              <description>
                This option will automatically determine the mininum/maximum
                This option will automatically determine the minimum/maximum
                latitude/longitude from the input camera model cube specified
                using the FROM parameter.
              </description>
@@ -538,7 +545,7 @@
           <option value="MAP">
              <brief> Read default range from map file</brief>
              <description>
                This option will read the mininum/maximum latitude/longitude
                This option will read the minimum/maximum latitude/longitude
                from the input map file.  All four values are expected to be
                defined.
              </description>
@@ -552,7 +559,7 @@
        <internalDefault>Use default range</internalDefault>
        <description>
          The minimum latitude of the ground range.   If this is entered by the
          user it will override the default camera or map file value.  By default,
          user, it will override the default camera or map file value.  By default,
          planetocentric latitudes are assumed unless the map file specifies
          otherwise.
        </description>
@@ -567,7 +574,7 @@
        <internalDefault>Use default range</internalDefault>
        <description>
          The maximum latitude of the ground range.   If this is entered by the
          user it will override the default camera or map file value.  By default,
          user, it will override the default camera or map file value.  By default,
          planetocentric latitudes are assumed unless the map file specifies
          otherwise.
        </description>
@@ -582,7 +589,7 @@
        <internalDefault>Use default range</internalDefault>
        <description>
          The minimum longitude of the ground range.   If this is entered by
          the user it will override the default camera or map value.  By
          the user, it will override the default camera or map value.  By
          default, positive east longitudes in the range of 0 to 360 are
          assumed unless the map file specifies otherwise.
        </description>
@@ -594,7 +601,7 @@
        <internalDefault>Use default range</internalDefault>
        <description>
          The maximum longitude of the ground range.   If this is entered by
          the user it will override the default camera or map value.  By
          the user, it will override the default camera or map value.  By
          default, positive east longitudes in the range of 0 to 360 are
          assumed unless the map file specifies otherwise.
        </description>
+10 −0
Original line number Diff line number Diff line
APPNAME = nocam2map

include $(ISISROOT)/make/isismake.tsts

commands:
	$(APPNAME) from=$(INPUT)/trdr9x9.cub \
	latcub=$(INPUT)/ddr9x9.cub+1 \
	loncub=$(INPUT)/ddr9x9.cub+2 \
	to=$(OUTPUT)/nocam2maptruth.cub > /dev/null;