Commit 9a0e5009 authored by Stuart Sides's avatar Stuart Sides
Browse files

Adding new gap capability for Cahandrayann1M3 and formatting for some base files

git-svn-id: http://subversion.wr.usgs.gov/repos/prog/isis3/branches/ipce@7222 41f8697f-d340-4b68-9986-7bafba869bb8
parent 8023db9c
Loading
Loading
Loading
Loading
+13 −0
Original line number Diff line number Diff line
@@ -42,6 +42,7 @@ namespace Isis {
    p_band(0), p_nbands(0), p_npixels(0), p_buf(0),
    p_pixelType(None), p_rawbuf(0) { }


  /**
   * Creates a Buffer object.  Note that the number of pixels in the buffer
   * can not exceed 2GB.
@@ -78,6 +79,7 @@ namespace Isis {
    Allocate();
  }


  //! Destroys the Buffer object and frees shape buffer.
  Buffer::~Buffer() {
    try {
@@ -93,6 +95,7 @@ namespace Isis {
    }
  }


  /**
   * @brief Assign the entire buffer to a constant double value
   * @param d Value to assign to the buffer
@@ -106,6 +109,7 @@ namespace Isis {
    return (*this);
  }


  /**
   * This method is used to set the base position of the shape buffer.  In general
   * it is used by BufferManager objects to progress sequentially through
@@ -140,6 +144,7 @@ namespace Isis {
    return (index % p_nsamps) + p_sample;
  }


  /**
   * Returns the line position associated with a shape buffer index.  The
   * shape buffer is one dimensional.  Let us assume a nsamps=2, nlines=3, and
@@ -158,6 +163,7 @@ namespace Isis {
    return sub_index / p_nsamps + p_line;
  }


  /**
   * Returns the band position associated with a shape buffer index.  The
   * shape buffer is one dimensional.  Let us assume a nsamps=2, nlines=3, and
@@ -173,6 +179,7 @@ namespace Isis {
    return index / (p_nsamps * p_nlines) + p_band;
  }


  /**
   * Returns the sample, line, and band position associated with a
   * shape buffer index.  Performs the same function as the Sample, Line,
@@ -190,6 +197,7 @@ namespace Isis {
    i_band = Band(index);
  }


  /**
   * Given a sample, line, and band position, this returns the appropriate index
   * in the shape buffer.
@@ -226,6 +234,7 @@ namespace Isis {
    return (index);
  }


  /**
   * Returns the value in the shape buffer at the given index.
   *
@@ -248,6 +257,7 @@ namespace Isis {
    return p_buf[index];
  }


  /**
   * Allows copying of the buffer contents to another Buffer
   *
@@ -278,6 +288,7 @@ namespace Isis {
    }
  }


  /**
   * Allows copying of the buffer contents of a larger buffer to another same size or smaller
   *   Buffer, using their base positions to relate data. This does not copy the raw buffer.
@@ -313,6 +324,7 @@ namespace Isis {
    return isSubareaOfIn;
  }


  /**
   * The copy constructor. Allows a new Buffer object to be created using
   * an existing Buffer object.
@@ -333,6 +345,7 @@ namespace Isis {
    Copy(rhs);
  }


  /**
   * Size or resize the memory buffer.
   *
+5 −3
Original line number Diff line number Diff line
@@ -152,6 +152,7 @@ namespace Isis {
               (BigInt)((p_maxBands - 1) / p_binc + 1));
  }


  /**
   * Sets the offset of the buffer. By default (if this method is not invoked)
   * the offsets are (0,0,0). Offsets are applied when computing the top-left
@@ -179,6 +180,7 @@ namespace Isis {
    p_boff = boff;
  }


  /**
   * Sets the position of the shape in the cube.  This shape fits
   * into the cube a specific number of times.
+76 −72
Original line number Diff line number Diff line
@@ -189,9 +189,11 @@ namespace Isis {
      *  @returns bool  true if the solution converged, false otherwise
      */
    template <typename Functor> 
    static bool brentsRootFinder(Functor &func, const QList<double> pt1, 
                          const QList<double> pt2, double tol, int maxIter, double &root) {
      double a=pt1[0], b=pt2[0], c, d=0, fa = pt1[1], fb = pt2[1], fc, p, q, r, s, t,tol1, bnew, fbnew, temp, deltaI,deltaF;
    static bool brentsRootFinder(Functor &func, const QList<double> pt1, const QList<double> pt2, 
                                 double tol, int maxIter, double &root) {
      double a = pt1[0], b = pt2[0], c, d = 0;
      double fa = pt1[1], fb = pt2[1], fc;
      double p, q, r, s, t, tol1, bnew, fbnew, temp, deltaI, deltaF;

      // offset used for improved numerical stability
      double offset = (a + b) / 2.0;
@@ -285,7 +287,7 @@ namespace Isis {
          throw IException(e, IException::Programmer, msg, _FILEINFO_);
        }
        d = c;  // thus d always equals the best guess from two iterations ago
        c = b;  //thus c always equals the best giess from the previous iteration
        c = b;  // thus c always equals the best guess from the previous iteration
        fc = fb;
        
        if ( (fa > 0) - (fa < 0) == (fbnew > 0) - (fbnew < 0) ) {
@@ -340,16 +342,18 @@ namespace Isis {
    private:
      /** Constructor
       *
       * This is private and to left undefined so this class cannot be instaniated
       * This is private and left undefined so this class cannot be instaniated. The functions 
       * are static and are intended to be used without an instance of the class. 
       */
      FunctionTools(); //no definition to be supplied
      FunctionTools();
     
      /** destructor
       *
       * This is private and to left undefined so this class cannot be instaniated
       */
      ~FunctionTools(); //no definition to be supplied
   }; //end FuntionTools class definition
      ~FunctionTools();

  }; // End FuntionTools class definition

}; // End namespace Isis

+67 −62
Original line number Diff line number Diff line
@@ -62,7 +62,7 @@ class LineOffsetFunctor :

    LineOffsetFunctor(Isis::Camera *camera, const Isis::SurfacePoint &surPt) {
      m_camera = camera;
      surfacePoint = surPt;
      m_surfacePoint = surPt;
    }


@@ -96,7 +96,7 @@ class LineOffsetFunctor :
      m_camera->Sensor::setTime(et);
 
      // Set ground
      if (!m_camera->Sensor::SetGround(surfacePoint, false)) {
      if (!m_camera->Sensor::SetGround(m_surfacePoint, false)) {
        IString msg = "Sensor::SetGround failed for surface point in LineScanCameraGroundMap.cpp"
                      " LineOffsetFunctor";
        throw IException(IException::Programmer, msg, _FILEINFO_);
@@ -127,7 +127,7 @@ class LineOffsetFunctor :

      // Try to use SetUndistortedFocalPlane, if that does not work use the distorted x,y
      // under the assumption (bad|good) that extrapolating the distortion
      // is causing the distorted x to be way off the sensor, and thus not very good anyway.
      // is causing the distorted x,y to be way off the sensor, and thus not very good anyway.
      if (m_camera->DistortionMap()->SetUndistortedFocalPlane(ux, uy)) {
        // Get the natural (distorted focal plane coordinates)
        dx = m_camera->DistortionMap()->FocalPlaneX();
@@ -151,7 +151,7 @@ class LineOffsetFunctor :


  private:
    SurfacePoint surfacePoint;
    SurfacePoint m_surfacePoint;
    Camera* m_camera;
};

@@ -306,17 +306,18 @@ namespace Isis {
    LineOffsetFunctor offsetFunc(p_camera,surfacePoint);
    SensorSurfacePointDistanceFunctor distanceFunc(p_camera,surfacePoint);

    /*********************************************************************************************
    if an approximate point is given use that as a start point for the secant method root search
    *********************************************************************************************/
    if (approxLine > 0) {  
    // METHOD #1
    // Use the line given as a start point for the secant method root search. 
    if (approxLine >= 0.5) {

      // convert the approxLine to an approximate time and offset
      p_camera->DetectorMap()->SetParent(p_camera->ParentSamples() / 2, approxLine);
      p_camera->DetectorMap()->SetParent(p_camera->ParentSamples() / 2.0, approxLine);
      approxTime = p_camera->time().Et();
  
      approxOffset = offsetFunc(approxTime);

      if (fabs(approxOffset) < 1e-2) { //no need to iteratively improve this root, it's good enough
      // Check to see if there is no need to improve this root, it's good enough
      if (fabs(approxOffset) < 1e-2) { 
        p_camera->Sensor::setTime(approxTime);
        // check to make sure the point isn't behind the planet
        if (!p_camera->Sensor::SetGround(surfacePoint, true)) {
@@ -370,8 +371,7 @@ namespace Isis {
          fh = f;
        }

        //See if we converged on the point so set up the undistorted
        //  focal plane values and return
        // See if we converged on the point so set up the undistorted focal plane values and return
        if (fabs(f) < 1e-2) {
          p_camera->Sensor::setTime(approxTime);
          // check to make sure the point isn't behind the planet
@@ -387,14 +387,13 @@ namespace Isis {

          return Success;       
        }
      }
      return Failure;
    }
      } // End itteration using a guess
      // return Failure; // Removed to let the lagrange method try to find the line if secant fails
    } // End use a guess


    /**********************************************************************************************
    no estimate given for the approximate line--quadratic approximation root finding
    **********************************************************************************************/
    // METHOD #2
    // The guess or middle line did not work so try estimating with a quadratic
    // The offsets are typically quadratic, so three points will be used to approximate a quadratic
    // as a first order attempt to find the root location(s)

@@ -460,6 +459,8 @@ namespace Isis {
    // Now that we have the coefficients of the quadratic look for roots 
    // (see Numerical Recipes Third Edition page 227)
    temp = quadPoly[1] * quadPoly[1] - 4.0 * quadPoly[0] * quadPoly[2];  //discriminant

    // THIS IS A PREMATURE FAILURE RETURN. IT SHOULD TRY THE NEXT METHON BEFORE FAILING
    if (temp < 0.0) {
      return Failure;  // there are apparently not any real roots on this image
    }
@@ -491,12 +492,13 @@ namespace Isis {
      root[i] = root[i]/scale + timeAverage;
    }

    // THIS IS A PREMATURE FAILURE RETURN. IT SHOULD TRY THE NEXT METHON BEFORE FAILING
    if (root.size() == 0) {
      return Failure;  // there are apparently not any roots on this image
    }

    // At the time of this writing ISIS made no attempt to support any sensors that were not "1 to 1".
    // Meaning that imaged the same point on the ground in multiple lines of the image
    // Meaning they imaged the same point on the ground in multiple lines of the image
    // therefore we must somehow reduce multiple possible roots to a single one,  the legacy 
    // code (replaced with this code) did this based on distance from the sensor to the target
    // the shortest distance being the winner.  For legacy consistency I have used the same logic below.  
@@ -507,13 +509,15 @@ namespace Isis {
    }

    // Save the root with the smallest dist
    {
      int j=0;   
    for (int i=1, j=0; i<root.size(); i++) {
      for (int i=1; i<root.size(); i++) {
        if (dist[i] < dist[j]) j=i;
      }      
        
      approxTime = root[j];  // Now we have our time
      approxOffset = offset[j];  // The offsets are saved to avoid recalculating it later
    }

    if (fabs(approxOffset) < 1.0e-2) { // No need to iteratively improve this root, it's good enough
      p_camera->Sensor::setTime(approxTime);
@@ -533,10 +537,9 @@ namespace Isis {
      return Success;
    }

    /**********************************************************************************************
    no estimate given for the approximate line, quadratic approximation insufficient, use Brent's
    method
    **********************************************************************************************/

    // METHOD #3
    // Estimated line and quadratic approximation insufficient, try Brent's method
    // The offsets are typically quadratic, so three points will be used to approximate a quadratic
    // as a first order attempt to find the root location(s)

@@ -577,7 +580,6 @@ namespace Isis {
        if (FunctionTools::brentsRootFinder <LineOffsetFunctor> (offsetFunc, pts[i-1], pts[i],
                                                                 1.0e-3, 200, temp)) {
          root << temp;

        }
      }      
    }
@@ -605,11 +607,14 @@ namespace Isis {
    }

    // Save the root with the smallest dist
    for (int i=1, j=0; i<root.size(); i++) {
    {
      int j=0;   
      for (int i=1; i<root.size(); i++) {
        if (dist[i] < dist[j]) j=i;
      }

      p_camera->Sensor::setTime(root[j]);
    }

    // No need to make sure the point isn't behind the planet, it was done above
    p_camera->Sensor::LookDirection(lookC);
+94 −26
Original line number Diff line number Diff line
@@ -65,6 +65,7 @@ void importImage(QString outputParamName, ProcessImportPds::PdsFileType fileType
  g_oBuff = NULL;
  g_totalLinesAdded = 0;
  g_utcTable = NULL;
  double calcOutputLines = 0;

  ProcessImportPds importPds;
  importPds.Progress()->SetText((QString)"Writing " + outputParamName + " file");
@@ -100,10 +101,13 @@ void importImage(QString outputParamName, ProcessImportPds::PdsFileType fileType
    // acquisition processing needs to be removed. There are four possible flip/mirror mode 
    // combinations.
    // 1.  Descending yaw / Forward orbit limb - No changes in sample or line
    // 2.  Descending yaw / Reverse orbit limb - Samples are reversed, first sample on west side of image
    // 3.  Ascending yaw / Forward orbit limb - Lines/times are reversed so northernmost image line first,
    //                                          Samples are reversed, first sample on west side of image
    // 4.  Ascending yaw / Reverse orbit limb - Lines/times are reversed so northernmost image line first,
    // 2.  Descending yaw / Reverse orbit limb - Samples are reversed, first sample on west side 
    //                                           of image
    // 3.  Ascending yaw / Forward orbit limb - Lines/times are reversed so northernmost image 
    //                                          line first, Samples are reversed, first sample on
    //                                          west side of image
    // 4.  Ascending yaw / Reverse orbit limb - Lines/times are reversed so northernmost image 
    //                                          line first,
    QString yawDirection = (QString) pdsLabel["CH1:SPACECRAFT_YAW_DIRECTION"];
    QString limbDirection = (QString) pdsLabel["CH1:ORBIT_LIMB_DIRECTION"];
    samplesNeedFlipped = ( ((yawDirection == "REVERSE") && (limbDirection == "DESCENDING")) ||
@@ -128,16 +132,37 @@ void importImage(QString outputParamName, ProcessImportPds::PdsFileType fileType

        QString instMode = (QString) pdsLabel["INSTRUMENT_MODE_ID"];
        // Initialize to the value for a GLOBAL mode observation
        g_expectedLineRate = .10176;
        g_expectedLineRate = 0.10176;
        if (instMode == "TARGET") {
          g_expectedLineRate = .05088;
          g_expectedLineRate = 0.05088;
        }

        iTime firstEt((QString)(*g_utcTable)[0]["UtcTime"]);
        iTime lastEt((QString)(*g_utcTable)[importPds.Lines()-1]["UtcTime"]);

        // The UTC line time table has been flipped in the same manner as the image lines, thus fabs
        outputLines = ceil(fabs(lastEt - firstEt) / g_expectedLineRate + 1.0);
        // Search the time table for gaps to come up with an output cube number of lines
        // The times in the table are documented as the time at the center of the exposure/frame, so
        // consecutive records in the time table should differ by the exposure rate, if not then
        // there is a potential gap.
        // This was calculated in the previous version of this code, but there is a minor difference 
        // between the calculation and the following brute force method.
        outputLines = 0;
        for (int rec = 0; rec < g_utcTable->Records() - 1; rec++) {
          outputLines++; // One for this line

          iTime thisEt((QString)(*g_utcTable)[rec]["UtcTime"]); 
          iTime nextEt((QString)(*g_utcTable)[rec+1]["UtcTime"]);
          double delta = fabs(nextEt - thisEt); // Time table may be assending or decenting times

          while (delta > g_expectedLineRate * 1.9) {
            outputLines++; // Big enough gap to need more line(s)
            delta -= g_expectedLineRate;
          }
        }
        outputLines++; // One more for the last line

        iTime firstEt((QString)(*g_utcTable)[0]["UtcTime"]); 
        iTime lastEt((QString)(*g_utcTable)[g_utcTable->Records()-1]["UtcTime"]);
        calcOutputLines = fabs((lastEt + g_expectedLineRate / 2.0) - 
                               (firstEt - g_expectedLineRate / 2.0)) / g_expectedLineRate;
      }
      else {
        QString msg = "Input file [" + in.expanded() +
@@ -147,6 +172,7 @@ void importImage(QString outputParamName, ProcessImportPds::PdsFileType fileType
    }
    else {
      outputLines = importPds.Lines();
      calcOutputLines = outputLines;
    }

    // Since the output cube possibly has more lines then the input PDS image, due to dropped 
@@ -168,14 +194,18 @@ void importImage(QString outputParamName, ProcessImportPds::PdsFileType fileType
    }
    else {
      importPds.StartProcess(writeCubeWithDroppedLines);
      g_results += PvlKeyword("LinesAdded", toString(g_totalLinesAdded));
      g_results += PvlKeyword("LinesFlipped", toString(linesNeedFlipped));
      g_results += PvlKeyword("SamplesFlipped", toString(samplesNeedFlipped));
      g_results += PvlKeyword("LinesAdded", toString(g_totalLinesAdded));
      g_results += PvlKeyword("OutputLines", toString(outputLines));
      g_results += PvlKeyword("CalculatedOutputLines", toString(calcOutputLines));
    }

    delete g_oBuff;

    // If the image lines need flipped then so does the UTC table, if it exists.
    // This does not need to be done before the main processing because the flipping of 
    // the image is done after the main processing.
    if (fileType != ProcessImportPds::L0) {
      if (linesNeedFlipped) {
        flipUtcTable(*g_utcTable);
@@ -213,6 +243,8 @@ void importImage(QString outputParamName, ProcessImportPds::PdsFileType fileType
}


// Processing function for writing all input PDS lines to the output cube.
// No dropped lines are inserted.
void writeCube(Buffer &in) {

  for (int i = 0; i < in.size(); i++) {
@@ -224,7 +256,8 @@ void writeCube(Buffer &in) {
}



// Processing function for writing all input PDS lines to the output cube with dropped lines
// inserted where the time table shows gaps.
void writeCubeWithDroppedLines(Buffer &in) {

  // Always write the current line to the output cube first
@@ -236,12 +269,16 @@ void writeCubeWithDroppedLines(Buffer &in) {
  (*g_oBuff)++;

  // Now check the UTC_TIME table and see if there is a gap (missing lines) after the TIME record
  // for the current line. If there are, add as many line as are necessary to fill the gap. 
  // for the current line. If there are, add as many lines as are necessary to fill the gap. 
  // Since the PDS files are in BIL order we are writeing to the ISIS cube in that order, so we 
  // do not need to write NULL lines to fill a gap until we have written the last band of line N,
  // and we don't have to check for gaps after the last lines of the PDS file. 

  if (in.Band() == g_oCube->bandCount() && in.Line() < g_utcTable->Records()) {

    QString tt = (QString)(*g_utcTable)[in.Line() - 1]["UtcTime"];
    QString ttt = (QString)(*g_utcTable)[in.Line()]["UtcTime"];

    iTime thisEt((QString)(*g_utcTable)[in.Line() - 1]["UtcTime"]);
    iTime nextEt((QString)(*g_utcTable)[in.Line()]["UtcTime"]);

@@ -270,8 +307,7 @@ void writeCubeWithDroppedLines(Buffer &in) {
}




// Transfere the needed PDS labels to the ISIS Cube and update them where necessary
void translateChandrayaan1M3Labels(Pvl& pdsLabel, Cube *ocube, Table& utcTable,
                                   ProcessImportPds::PdsFileType fileType) {
  Pvl outLabel;
@@ -313,21 +349,53 @@ void translateChandrayaan1M3Labels(Pvl& pdsLabel, Cube *ocube, Table& utcTable,

    SpiceInt sclkCode = -86;

    QString startTime = inst["SpacecraftClockStartCount"];
    double et;
    scs2e_c(sclkCode, startTime.toLatin1().data(), &et);
    iTime startEt(et);
    inst.findKeyword("StartTime").setValue(startEt.UTC());
    // Remmoved when we found out the lable counts are not as correct as we need. We use the time
    // tables instead (see below)
    //QString startTime = inst["SpacecraftClockStartCount"];
    //double et;
    //scs2e_c(sclkCode, startTime.toAscii().data(), &et);
    //iTime startEt(et);
    //inst.findKeyword("StartTime").setValue(startEt.UTC());

    //QString stopTime = inst["SpacecraftClockStopCount"];
    //scs2e_c(sclkCode, stopTime.toAscii().data(), &et);
    //iTime stopEt(et);
    //inst.findKeyword("StopTime").setValue(stopEt.UTC());

    // Replace code above with this
    // The start and stop times in the PDS labels do not match the UTC table times.
    // Assume the UTC table times are better, so change the labels to match the table
    // The start and stop clock counts need to match the start/stop time, so convert the times
    // to new clock counts.
    iTime firstEt((QString)(*g_utcTable)[0]["UtcTime"]);
    iTime lastEt((QString)(*g_utcTable)[utcTable.Records()-1]["UtcTime"]);

    QString stopTime = inst["SpacecraftClockStopCount"];
    scs2e_c(sclkCode, stopTime.toLatin1().data(), &et);
    iTime stopEt(et);
    inst.findKeyword("StopTime").setValue(stopEt.UTC());
    // The table is in assending order
    // The table contains the middle of the exposure. include times to cover the beginning of
    // line 1 and the end of line NL
    if (firstEt < lastEt) {
      firstEt = firstEt - (g_expectedLineRate / 2.0);
      lastEt = lastEt + (g_expectedLineRate / 2.0);
    }
    else {
      firstEt = lastEt - (g_expectedLineRate / 2.0);
      lastEt = firstEt + (g_expectedLineRate / 2.0);
    }

  ocube->putGroup(inst);
    inst.findKeyword("StartTime").setValue(firstEt.UTC());
    SpiceChar startClockString[100];
    sce2s_c (sclkCode, firstEt.Et(), 100, startClockString);
    QString startClock(startClockString);
    inst.findKeyword("SpacecraftClockStartCount").setValue(startClock);

    inst.findKeyword("StopTime").setValue(lastEt.UTC());
    SpiceChar stopClockString[100];
    sce2s_c (sclkCode, lastEt.Et(), 100, stopClockString);
    QString stopClock(stopClockString);
    inst.findKeyword("SpacecraftClockStopCount").setValue(stopClock);
  }

  ocube->putGroup(inst);

  if (fileType == ProcessImportPds::L0 || fileType == ProcessImportPds::Rdn) {
    // Setup the band bin group
Loading