Commit d2a31f0c authored by Stuart Sides's avatar Stuart Sides
Browse files

First commit of Chandrayaan M3 code. Ref#1553, Ref#167701678

git-svn-id: http://subversion.wr.usgs.gov/repos/prog/isis3/trunk@5761 41f8697f-d340-4b68-9986-7bafba869bb8
parent cf9cbf3b
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -73,7 +73,7 @@ plugins:


clean:
	rm -f lib/lib*.so* lib/lib*.dylib
	rm -f lib/lib*.so* lib/lib*.dylib lib/lib*.a
	cd lib && rm -rf *.framework
	@for plugs in plugins/*; do \
	  if [ -d $$plugs -a $$plugs != "plugins/CVS" ]; \
+1 −2
Original line number Diff line number Diff line
@@ -78,7 +78,7 @@ EndGroup
########################################################

Group = SessionLog
  TerminalOutput = On
  TerminalOutput = Off
  FileOutput     = On
  FileName       = print.prt
  FileAccess     = Append
@@ -158,7 +158,6 @@ Group = DataDirectory
  Apollo17     = $ISIS3DATA/apollo17
  Base         = $ISIS3DATA/base
  Cassini      = $ISIS3DATA/cassini
  Chan1        = $ISIS3DATA/chan1
  Chandrayaan1 = $ISIS3DATA/chandrayaan1
  Clementine1  = $ISIS3DATA/clementine1
  Control      = $ISIS3DATA/control
+4 −0
Original line number Diff line number Diff line
@@ -22,6 +22,8 @@
 */
#include "CameraPointInfo.h"

#include <iomanip>

#include "Brick.h"
#include "Camera.h"
#include "Cube.h"
@@ -90,6 +92,8 @@ namespace Isis {
                                      const bool outside, const bool errors) {
    if (CheckCube()) {
      bool passed = m_camera->SetImage(sample, line);


      return GetPointInfo(passed, outside, errors);
    }
    // Should never get here, error will be thrown in CheckCube()
+25 −8
Original line number Diff line number Diff line
@@ -401,14 +401,31 @@ namespace Isis {
                    + m_tableName +"]. The PDS label file is probably invalid";
      throw IException(IException::Unknown, msg.toStdString(), _FILEINFO_);
    }
    m_recordBytes = (int) label.findKeyword("RECORD_BYTES");
    //  Get some pertinent information from the label
    PvlObject &tabobj = label.findObject(m_tableName);
    PvlObject &tabObj = label.findObject(m_tableName);
    if (tabObj.hasKeyword("RECORD_BYTES")) {
      m_recordBytes = (int) tabObj.findKeyword("RECORD_BYTES");
    }
    else {
      m_recordBytes = (int) label.findKeyword("RECORD_BYTES");
    }
    
    QString trueTableName;
    PvlObject *tableDetails = &tabObj;
    if (label.hasKeyword("^" + m_tableName)) {
      trueTableName = m_tableName;
      pdsTableFile = FileName(pdsLabFile).path() + "/" 
                     + label["^" + m_tableName][0];
    m_trows = (int) tabobj.findKeyword("ROWS");
    int ncols =  (int) tabobj.findKeyword("COLUMNS");
    m_pdsTableType = QString(tabobj.findKeyword("INTERCHANGE_FORMAT"));
    }
    else if (tabObj.objects() == 1) {
      trueTableName = tabObj.object(0).name();
      tableDetails = &tabObj.object(0);
      pdsTableFile = FileName(pdsLabFile).path() + "/" 
                     + tabObj["^" + trueTableName][0];
    }
    m_trows = (int) tableDetails->findKeyword("ROWS");
    int ncols =  (int) tableDetails->findKeyword("COLUMNS");
    m_pdsTableType = QString(tableDetails->findKeyword("INTERCHANGE_FORMAT"));
    if (m_pdsTableType != "ASCII" && m_pdsTableType.toUpper() != "BINARY") {
      QString msg = "Unable to import the PDS table [" + m_tableName 
                    + "] from the PDS file [" 
@@ -417,12 +434,12 @@ namespace Isis {
                    + "] is not supported. Valid values are ASCII or BINARY.";
      throw IException(IException::User, msg.toStdString(), _FILEINFO_);
    }
    m_rowBytes = tabobj.findKeyword("ROW_BYTES");
    m_rowBytes = tableDetails->findKeyword("ROW_BYTES");

    m_coldesc.clear();
    PvlObject::PvlObjectIterator colobj = tabobj.beginObject();
    PvlObject::PvlObjectIterator colobj = tableDetails->beginObject();
    int icol(0);
    while (colobj != tabobj.endObject()) {
    while (colobj != tableDetails->endObject()) {
      if (colobj->isNamed("COLUMN")) {
        m_coldesc.push_back(getColumnDescription(*colobj, icol));
        icol++;
+158 −64
Original line number Diff line number Diff line
@@ -65,15 +65,26 @@ class LineOffsetFunctor :
    }
    ~LineOffsetFunctor() {}
            
         
       
    /** Compute the number of lines between the current line (i.e., the line imaged at the et as set
     *  in the camera model) and the line number where the argument et would hit the focal
     *  plane.
     *  
     * @param et The et at the new postion
     *
     * @return Line off (see description)
     */
    double operator()(double et) {
      double lookC[3] = {0.0, 0.0, 0.0};
      double ux = 0.0, uy = 0.0;
      double dx = 0.0, dy = 0.0;
      double ux = 0.0;
      double uy = 0.0;
      double dx = 0.0;
      double dy = 0.0;

      //verify the time is with the cache bounds
      // Verify the time is with the cache bounds
      double startTime = m_camera->cacheStartTime().Et();
      double endTime = m_camera->cacheEndTime().Et();
      //std::pair <double, double> cacheTimeBounds = m_camera->StartEndEphemerisTimes();
      if (et < startTime || et > endTime) {
        IString msg = "Ephemeris time passed to LineOffsetFunctor is not within the image "
                      "cache bounds";
@@ -82,28 +93,48 @@ class LineOffsetFunctor :
     
      m_camera->Sensor::setTime(et);
 
      //set ground
      // Set ground
      if (!m_camera->Sensor::SetGround(surfacePoint, false)) {
        IString msg = "Sensor::SetGround failed for surface point in LineScanCameraGroundMap.cpp"
                      " LineOffsetFunctor";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
   
      //calculate the undistorted focal plane coordinates
      // Calculate the undistorted focal plane coordinates
      m_camera->Sensor::LookDirection(lookC);
      ux = m_camera->FocalLength() * lookC[0] / lookC[2];
      uy = m_camera->FocalLength() * lookC[1] / lookC[2];

      //set the undistorted focal plane coordinates
      if (!m_camera->DistortionMap()->SetUndistortedFocalPlane(ux, uy)) {
        IString msg = "DistortionMap::SetUndistoredFocalPlane failed for surface point in "
                      "LineScanCameraGroundMap.cpp LineOffsetFunctor";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }

      //get the natural (distorted focal plane coordinates)
      // This was replaced with the code below to get Chandrayaan M3 to work.
      // SetUndistortedFocalPlane was failing a majority of the time, causing most SetGround call
      // to fail. Even when it did succeed, it was producing non-continous return values.
        // Set the undistorted focal plane coordinates
//        if (!m_camera->DistortionMap()->SetUndistortedFocalPlane(ux, uy)) {
//          IString msg = "DistortionMap::SetUndistoredFocalPlane failed for surface point in "
//                        "LineScanCameraGroundMap.cpp LineOffsetFunctor";
//          throw IException(IException::Programmer, msg, _FILEINFO_);
//        }

        // Get the natural (distorted focal plane coordinates)
//        dx = m_camera->DistortionMap()->FocalPlaneX();
//        dy = m_camera->DistortionMap()->FocalPlaneY();
//        std::cout << "use dist" << std::endl;
//      }


      // 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.
      if (m_camera->DistortionMap()->SetUndistortedFocalPlane(ux, uy)) {
        // Get the natural (distorted focal plane coordinates)
        dx = m_camera->DistortionMap()->FocalPlaneX();
        dy = m_camera->DistortionMap()->FocalPlaneY();
      }
      else {
        dx = ux;
        dy = uy;
      }

      if (!m_camera->FocalPlaneMap()->SetFocalPlane(dx, dy)) {
        IString msg = "FocalPlaneMap::SetFocalPlane failed for surface point in "
@@ -111,11 +142,12 @@ class LineOffsetFunctor :
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }     

      //return the offset
      // Return the offset
      return (m_camera->FocalPlaneMap()->DetectorLineOffset() - 
              m_camera->FocalPlaneMap()->DetectorLine());///m_camera->PixelPitch();
              m_camera->FocalPlaneMap()->DetectorLine());
    }


  private:
    SurfacePoint surfacePoint;
    Camera* m_camera;
@@ -136,12 +168,14 @@ class SensorSurfacePointDistanceFunctor :
    }
    ~SensorSurfacePointDistanceFunctor() {}
            
                
    double operator()(double et) {
      double s[3], p[3];

      //verify the time is with the cache bounds
      double startTime = m_camera->cacheStartTime().Et();
      double endTime = m_camera->cacheEndTime().Et();
      //std::pair <double, double> cacheTimeBounds = m_camera->StartEndEphemerisTimes();

      if (et < startTime || et > endTime) {
        IString msg = "Ephemeris time passed to SensorSurfacePointDistanceFunctor is not within the image "
                      "cache bounds";
@@ -166,6 +200,19 @@ class SensorSurfacePointDistanceFunctor :


namespace Isis {

  /** Constructor
   *
   * @param cam pointer to camera model
   */
  LineScanCameraGroundMap::LineScanCameraGroundMap(Camera *cam) : CameraGroundMap(cam) {}


  /** Destructor
   *
   */
  LineScanCameraGroundMap::~LineScanCameraGroundMap() {}

  /** Compute undistorted focal plane coordinate from ground position
   *
   * @param lat planetocentric latitude in degrees
@@ -176,13 +223,16 @@ namespace Isis {
  bool LineScanCameraGroundMap::SetGround(const Latitude &lat,
      const Longitude &lon) {
    Distance radius(p_camera->LocalRadius(lat, lon));

    if (radius.isValid()) {
      return SetGround(SurfacePoint(lat, lon, radius));
    } else {
    }
    else {
      return false;
    }
  }


  /** Compute undistorted focal plane coordinate from ground position
   *
   * @param lat planetocentric latitude in degrees
@@ -213,17 +263,20 @@ namespace Isis {
    return false;
  }


  double LineScanCameraGroundMap::FindSpacecraftDistance(int line,
      const SurfacePoint &surfacePoint) {

    CameraDetectorMap *detectorMap = p_camera->DetectorMap();
    detectorMap->SetParent(p_camera->ParentSamples() / 2, line);
    if(!p_camera->Sensor::SetGround(surfacePoint, false)) return DBL_MAX;
    if (!p_camera->Sensor::SetGround(surfacePoint, false)) {
      return DBL_MAX;
    }

    return p_camera->SlantDistance();
  }



  LineScanCameraGroundMap::FindFocalPlaneStatus
      LineScanCameraGroundMap::FindFocalPlane(const int &approxLine,
                                              const SurfacePoint &surfacePoint) {
@@ -231,14 +284,15 @@ namespace Isis {
    //CameraDistortionMap *distortionMap = p_camera->DistortionMap();
    //CameraFocalPlaneMap *focalMap = p_camera->FocalPlaneMap();

    double approxTime=0,approxOffset=0;
    double approxTime=0;
    double approxOffset=0;
    double lookC[3] = {0.0, 0.0, 0.0};
    double ux = 0.0, uy = 0.0;
    double ux = 0.0;
    double uy = 0.0;
    //double dx = 0.0, dy = 0.0;
    //double s[3], p[3];
    const double cacheStart = p_camera->Spice::cacheStartTime().Et();
    const double cacheEnd = p_camera->Spice::cacheEndTime().Et();
    int j=0;   

    double lineRate = ((LineScanCameraDetectorMap *)p_camera->DetectorMap())->LineRate(); //line rate

@@ -277,8 +331,12 @@ namespace Isis {

      //starting times for the secant method, kept within the domain of the cache
      xh = approxTime;
      if (xh + lineRate < cacheEnd) xl = xh + lineRate;
        else                        xl = xh - lineRate;
      if (xh + lineRate < cacheEnd) {
        xl = xh + lineRate;
      }
      else {
        xl = xh - lineRate;
      }

      //starting offsets
      fh = approxOffset;  //the first is already calculated
@@ -286,7 +344,9 @@ namespace Isis {

      // Iterate to refine the given approximate time that the instrument imaged the ground point
      for (int j=0; j < 10; j++) {
        if ( fl-fh == 0.0) return Failure;
        if (fl-fh == 0.0) {
          return Failure;
        }
        double etGuess = xl + (xh - xl) * fl / (fl - fh);

        if (etGuess < cacheStart) etGuess = cacheStart;
@@ -330,12 +390,14 @@ namespace Isis {
    /**********************************************************************************************
    no estimate given for the approximate line--quadratic approximation root finding
    **********************************************************************************************/
    //the offsets are typically quadratic, so three points will be used to approximate 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)

    //the three nodes to be used to approximate the quadratic
    // The three nodes to be used to approximate the quadratic
    double offsetNodes[3];
    double timeNodes[3],timeAverage,scale;
    double timeNodes[3];
    double timeAverage;
    double scale;
    QList<double> root;
    QList<double> offset;
    QList<double> dist;
@@ -344,7 +406,8 @@ namespace Isis {
    timeNodes[2] = cacheEnd;
    timeNodes[1] = (cacheStart+cacheEnd) / 2.0; //middle time

    double quadPoly[3],temp;
    double quadPoly[3];
    double temp;

    for (int i=0; i<3; i++) {
      offsetNodes[i] = offsetFunc(timeNodes[i]);
@@ -356,7 +419,10 @@ namespace Isis {
    timeNodes[1] -= timeAverage;
    timeNodes[2] -= timeAverage;

    scale = 1.0/sqrt( (timeNodes[0] - timeNodes[2])*(timeNodes[0] - timeNodes[2]) + (offsetNodes[0] - offsetNodes[2])*(offsetNodes[0] - offsetNodes[2]) );
    scale = 1.0 / sqrt((timeNodes[0] - timeNodes[2]) * 
                       (timeNodes[0] - timeNodes[2]) + 
                       (offsetNodes[0] - offsetNodes[2]) * 
                       (offsetNodes[0] - offsetNodes[2]));

    timeNodes[0] *= scale;
    timeNodes[1] *= scale;
@@ -366,8 +432,9 @@ namespace Isis {
    offsetNodes[1] *= scale;
    offsetNodes[2] *= scale;

    //use lagrange interpolating polynomials to find the coefficients of the quadratic, 
    //  there are many ways to do this; I chose to do it this way because it is pretty straight forward and cheap
    // Use lagrange interpolating polynomials to find the coefficients of the quadratic, 
    // there are many ways to do this; I chose to do it this way because it is pretty straight 
    // forward and cheap
    quadPoly[0] = quadPoly[1] = quadPoly[2] = 0.0;

    temp = offsetNodes[0] / ((timeNodes[0] - timeNodes[1]) * (timeNodes[0] - timeNodes[2]));
@@ -385,50 +452,72 @@ namespace Isis {
    quadPoly[1] += temp * (-timeNodes[0] - timeNodes[1]);
    quadPoly[2] += temp * timeNodes[0] * timeNodes[1];

    //now that we have the coefficients of the quadratic look for roots 
    // 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
    if (temp < 0.0) return Failure;  //there are apparently not any real roots on this image
    if (temp < 0.0) {
      return Failure;  //there are apparently not any real roots on this image
    }

    if(quadPoly[1] >= 0.0) temp = -0.5*(quadPoly[1] + sqrt(temp));
    else                   temp = -0.5*(quadPoly[1] - sqrt(temp));
    if (quadPoly[0] != 0.0) root.push_back(temp/quadPoly[0]);
    if (quadPoly[2] != 0.0) root.push_back(quadPoly[2]/temp);
    if (quadPoly[1] >= 0.0) {
      temp = -0.5 * (quadPoly[1] + sqrt(temp));
    }
    else {
      temp = -0.5 * (quadPoly[1] - sqrt(temp));
    }

    if (quadPoly[0] != 0.0) {
      root.push_back(temp/quadPoly[0]);
    }

    if (quadPoly[2] != 0.0) {
      root.push_back(quadPoly[2]/temp);
    }

    //check to see if the roots are in the time interval of the cache
    for (int i=root.size()-1; i>=0; i--) {
      if ( root[i] < timeNodes[0] || root[i] > timeNodes[2] ) 
      if ( root[i] < timeNodes[0] || root[i] > timeNodes[2] ) {
        root.removeAt(i);
      }
    }
     
    //return the calculated roots to the original system
    for (int i=0;i<root.size();i++) root[i] = root[i]/scale + timeAverage;
    for (int i=0; i<root.size(); i++) {
      root[i] = root[i]/scale + timeAverage;
    }

    if (root.size()==0) return Failure;  //there are apparently not any roots on this image
    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
    //  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
    // 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
    // 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.  


    for (int i=0; i<root.size(); i++) {  //Offset/dist calculation loop
      dist << distanceFunc(root[i]);
      offset << offsetFunc(root[i]);
    }  //end Offset/dist calculation loop
    //save the root with the smallest dist
    }

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

    if (fabs(approxOffset) < 1e-2) { //no need to iteratively improve this root, it's good enough
    if (fabs(approxOffset) < 1.0e-2) { // No need to iteratively improve this root, it's good enough
      p_camera->Sensor::setTime(approxTime);
      //check to make sure the point isn't behind the planet

      // Check to make sure the point isn't behind the planet
      if (!p_camera->Sensor::SetGround(surfacePoint, true)) {
        return Failure;
      }
       
      p_camera->Sensor::LookDirection(lookC);
      ux = p_camera->FocalLength() * lookC[0] / lookC[2];
      uy = p_camera->FocalLength() * lookC[1] / lookC[2];
@@ -443,7 +532,7 @@ namespace Isis {
    no estimate given for the approximate line, quadratic approximation insufficient, use Brent's
    method
    **********************************************************************************************/
    //the offsets are typically quadratic, so three points will be used to approximate 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)

    // The above sections are sufficient for finding the correct times for the vast majority of 
@@ -456,8 +545,7 @@ namespace Isis {
    // that are bracketed by the five points defined in the quadratic solution method above.
    // The root with the shortest distance to the camera will be returned

    //Get everything ordered for iteration
    //  combine and sort the five already defined points
    // Get everything ordered for iteration combine and sort the five already defined points
    QList <QList <double> > pts;

    for (int i=0; i<3; i++) {
@@ -477,17 +565,19 @@ namespace Isis {
    qSort(pts.begin(), pts.end(), ptXLessThan);
    
    root.clear();
    
    for (int i=1; i<pts.size(); i++) {
      //if the signs of the two offsets are not the same they bracket at least one root
      // If the signs of the two offsets are not the same they bracket at least one root
      if ( (pts[i-1][1] > 0) - (pts[i-1][1] < 0) != (pts[i][1] > 0) - (pts[i][1] < 0) ) {
        double temp;
        if (FunctionTools::brentsRootFinder <LineOffsetFunctor> (offsetFunc,pts[i-1],pts[i],1.0e-3,200,temp))
        if (FunctionTools::brentsRootFinder <LineOffsetFunctor> (offsetFunc, pts[i-1], pts[i],
                                                                 1.0e-3, 200, temp)) {
          root << temp;

        }
      }      
    }

    //discard any roots that are looking through the planet
    // Discard any roots that are looking through the planet
    for (int i = root.size()-1; i>=0; i--) {
      p_camera->Sensor::setTime(root[i]);
      //check to make sure the point isn't behind the planet
@@ -496,24 +586,27 @@ namespace Isis {
      }
    }

    //if none of the roots remain...
    if (root.size() == 0) return Failure;

    // If none of the roots remain...
    if (root.size() == 0) {
      return Failure;
    }

    //choose from the remaining roots, the solution with the smallest distance to target
    // Choose from the remaining roots, the solution with the smallest distance to target
    dist.clear();
    offset.clear();
    for (int i=0; i<root.size(); i++) {  // Offset/dist calculation loop
      dist << distanceFunc(root[i]);
      offset << offsetFunc(root[i]);
    }  //end Offset/dist calculation loop
    //save the root with the smallest dist
    }

    // Save the root with the smallest dist
    for (int i=1, j=0; 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

    // No need to make sure the point isn't behind the planet, it was done above
    p_camera->Sensor::LookDirection(lookC);
    ux = p_camera->FocalLength() * lookC[0] / lookC[2];
    uy = p_camera->FocalLength() * lookC[1] / lookC[2];
@@ -525,6 +618,7 @@ namespace Isis {
  }
}


bool ptXLessThan(const QList<double> l1, const QList<double> l2) {
  return l1[0] < l2[0];
}
Loading