Commit 288dec17 authored by Jesse Mapel's avatar Jesse Mapel
Browse files

Added new app to check polynomial fit error, checkfit. 

git-svn-id: http://subversion.wr.usgs.gov/repos/prog/isis3/branches/PiecewisePolynomials@7870 41f8697f-d340-4b68-9986-7bafba869bb8
parent 7688bdb0

isis/src/base/apps/checkfit/Makefile

0 → 100644
+7 −0
Original line number Diff line number Diff line 
ifeq ($(ISISROOT), $(BLANK))

.SILENT:

error:

	echo "Please set ISISROOT";

else

	include $(ISISROOT)/make/isismake.apps

endif
 
 No newline at end of file

isis/src/base/apps/checkfit/checkfit.cpp

0 → 100644
+157 −0
Original line number Diff line number Diff line 
#include "Isis.h"



#include <cmath>

#include <iostream>



#include <QFile>

#include <QMap>

#include <QPair>

#include <QString>

#include <QTextStream>



#include "Camera.h"

#include "Cube.h"

#include "FileList.h"

#include "FileName.h"

#include "IException.h"

#include "PiecewisePolynomial.h"

#include "Progress.h"

#include "SpicePosition.h"

#include "SpiceRotation.h"

#include "UserInterface.h"



using namespace std;

using namespace Isis;



QPair<double, double> testFit(FileName inCubeFile, int positionDegree, int positionSegments,

                              int pointingDegree, int pointingSegments);



void IsisMain() {

  UserInterface &ui = Application::GetUserInterface();



  // Read in the list of cubes to check

  FileList cubeList;

  cubeList.read( ui.GetFileName("FROMLIST") );



  // Get the fit parameters

  int positionDegree = ui.GetInteger("SPKDEGREE");

  int positionSegments = ui.GetInteger("SPKSEGMENTS");

  int pointingDegree = ui.GetInteger("CKDEGREE");

  int pointingSegments = ui.GetInteger("CKSEGMENTS");



  // Setup the map for storing fit quality

  QMap<QString, QPair<double, double> > qualityMap;



  // Setup the progress tracker

  Progress cubeProgress;

  cubeProgress.SetMaximumSteps(cubeList.size());

  cubeProgress.CheckStatus();



  // Compute a test fit for each cube

  for (int cubeIndex = 0; cubeIndex < cubeList.size(); cubeIndex++) {

    FileName cubeFileName = cubeList[cubeIndex];

    QPair<double, double> fitQuality;

    try {

      cubeProgress.CheckStatus();

      fitQuality = testFit(cubeFileName,

                           positionDegree, positionSegments,

                           pointingDegree, pointingSegments);

    }

    catch(IException &e) {

      QString warning = "**WARNING** Failed checking cube [" + cubeFileName.expanded() + "].";

      std::cerr << warning << std::endl << e.toString() << std::endl;

      continue;

    }

    qualityMap.insert(cubeFileName.expanded(), fitQuality);

  }



  // Open the TO file for writing

  FileName outFileName = ui.GetFileName("TO");

  QFile outFile(outFileName.expanded());

  if (outFile.open(QFile::WriteOnly |QFile::Truncate)) {

    QTextStream outWriter(&outFile);

    // Output the header

    outWriter << "Cube, Position Error (km), Pointing Error (Rad)\n";

    QList<QString> cubeNames = qualityMap.keys();

    for (int i = 0; i < cubeNames.size(); i++) {

      QString cubeName = cubeNames[i];

      QPair<double, double> fitQuality = qualityMap.value(cubeName);

      outWriter << cubeName << ", "

                << toString(fitQuality.first) << ", "

                << toString(fitQuality.second) << "\n";

    }

  }

  else {

    QString msg = "Failed opening output file [" + outFileName.expanded() + "].";

    throw IException(IException::Io, msg, _FILEINFO_);

  }

}



QPair<double, double> testFit(FileName inCubeFile, int positionDegree, int positionSegments,

                              int pointingDegree, int pointingSegments) {

  Cube inCube(inCubeFile);

  Camera *inCam = inCube.camera();

  SpicePosition *instPosition = inCam->instrumentPosition();

  SpiceRotation *instRotation = inCam->instrumentRotation();



  // Fit the position

  instPosition->SetPolynomialDegree(positionDegree);

  instPosition->setPolynomialSegments(positionSegments);

  PiecewisePolynomial positionPoly = instPosition->testFit();



  // Fit the rotation

  instRotation->SetPolynomialDegree(pointingDegree);

  instRotation->setPolynomialSegments(pointingSegments);

  PiecewisePolynomial rotationPoly = instRotation->testFit();



  // Compute the position RMS

  double sumSquaredPositionError = 0;

  double positionBaseTime = instPosition->GetBaseTime();

  double positionTimeScale = instPosition->GetTimeScale();

  std::vector<double> positionSampleTimes = instPosition->timeCache();

  int positionSampleCount = positionSampleTimes.size();

  for (int i = 0; i < positionSampleCount; i++) {

    double error = 0;

    double scaledTime = (positionSampleTimes[i] - positionBaseTime) / positionTimeScale;

    std::vector<double> measuredCoord = instPosition->SetEphemerisTime(positionSampleTimes[i]);

    std::vector<double> estimatedCoord = positionPoly.evaluate(scaledTime);

    error += (measuredCoord[0] - estimatedCoord[0]) * (measuredCoord[0] - estimatedCoord[0]);

    error += (measuredCoord[1] - estimatedCoord[1]) * (measuredCoord[1] - estimatedCoord[1]);

    error += (measuredCoord[2] - estimatedCoord[2]) * (measuredCoord[2] - estimatedCoord[2]);

    sumSquaredPositionError += sqrt(error);

  }

  double positionRMS = sqrt(sumSquaredPositionError / positionSampleCount);



  // Compute the rotation RMS

  double sumSquaredRotationError = 0;

  double rotationBaseTime = instRotation->GetBaseTime();

  double rotationTimeScale = instRotation->GetTimeScale();

  std::vector<double> rotationSampleTimes = instRotation->timeCache();

  int rotationSampleCount = rotationSampleTimes.size();

  double start1 = 0.; // value of 1st angle1 in cache

  double start3 = 0.; // value of 1st angle3 in cache

  for (int i = 0; i < rotationSampleCount; i++) {

    double error = 0;

    double scaledTime = (rotationSampleTimes[i] - rotationBaseTime) / rotationTimeScale;

    instRotation->SetEphemerisTime(rotationSampleTimes[i]);

    std::vector<double> measuredAngles = instRotation->Angles(3, 1, 3);

    // Fix the angles crossing the domain bound

    if (i == 0) {

      start1 = measuredAngles[0];

      start3 = measuredAngles[2];

    }

    else {

      measuredAngles[0] = instRotation->WrapAngle(start1, measuredAngles[0]);

      measuredAngles[2] = instRotation->WrapAngle(start3, measuredAngles[2]);

    }

    std::vector<double> estimatedAngles = rotationPoly.evaluate(scaledTime);

    error += (measuredAngles[0] - estimatedAngles[0]) * (measuredAngles[0] - estimatedAngles[0]);

    error += (measuredAngles[1] - estimatedAngles[1]) * (measuredAngles[1] - estimatedAngles[1]);

    error += (measuredAngles[2] - estimatedAngles[2]) * (measuredAngles[2] - estimatedAngles[2]);

    sumSquaredRotationError += sqrt(error);

  }

  double rotationRMS = sqrt(sumSquaredRotationError / rotationSampleCount);



  return QPair<double, double>(positionRMS, rotationRMS);

}

isis/src/base/apps/checkfit/checkfit.xml

0 → 100644
+113 −0
Original line number Diff line number Diff line 
<?xml version="1.0" encoding="UTF-8"?>



<application name="checkfit" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:noNamespaceSchemaLocation="http://isis.astrogeology.usgs.gov/Schemas/Application/application.xsd">

  <brief>

    Test the quality of polynomial fits to cube exterior orientation

  </brief>



  <description>

    This program takes a list of cubes and tests the quality of the polynomial

    fit over instrument position and pointing values. This emulates how the

    jigsaw app converts exterior orientation into polynomials prior to

    adjustment. Poor quality fits can result in jigsaw failing.

  </description>



  <category>

    <categoryItem>Utility</categoryItem>

  </category>



  <seeAlso>

    <applications>

      <item>rotate</item>

      <item>flip</item>

    </applications>

  </seeAlso>



  <history>

    <change name="Jesse Mapel" date="2017-07-17">

      Original version

    </change>



  </history>



  <groups>

    <group name="Files">

      <parameter name="FROMLIST">

        <type>filename</type>

        <fileMode>input</fileMode>

        <brief>

          List of input cubes to check

        </brief>

        <description>

          The list of input cubes that a test polynomial fit will be done for.

        </description>

        <filter>

          *.lis

        </filter>

      </parameter>



      <parameter name="TO">

        <type>filename</type>

        <fileMode>output</fileMode>

        <brief>

          Output file that contains the quality of fits

        </brief>

        <description>

          This text file will contain the quality of fit for each input cube

          sorted from best fit to worst fit. The first column is the cube name,

          the second column is the quality of the instrument position fit, and

          the third column is the quality of the instrument pointing fit.

        </description>

        <filter>

          *.txt, *.csv

        </filter>

      </parameter>

    </group>



    <group name="Fit options">

      <parameter name="SPKDEGREE">

        <type>integer</type>

        <brief>

          The degree of the instrument position polynomial fit.

        </brief>

        <description>

          The degree of the instrument position polynomial fit.

        </description>

        <minimum inclusive="yes">0</minimum>

      </parameter>



      <parameter name="SPKSEGMENTS">

        <type>integer</type>

        <brief>

          The number of polynomial segments used in the instrument position fit.

        </brief>

        <description>

          The number of polynomial segments used in the instrument position fit.

        </description>

        <minimum inclusive="yes">1</minimum>

      </parameter>



      <parameter name="CKDEGREE">

        <type>integer</type>

        <brief>

          The degree of the instrument pointing polynomial fit.

        </brief>

        <description>

          The degree of the instrument pointing polynomial fit.

        </description>

        <minimum inclusive="yes">0</minimum>

      </parameter>



      <parameter name="CKSEGMENTS">

        <type>integer</type>

        <brief>

          The number of polynomial segments used in the instrument pointing fit.

        </brief>

        <description>

          The number of polynomial segments used in the instrument pointing fit.

        </description>

        <minimum inclusive="yes">1</minimum>

      </parameter>

    </group>

  </groups>

</application>

isis/src/base/objs/SpicePosition/SpicePosition.cpp

+97 −4
Original line number Diff line number Diff line
@@ -961,11 +961,29 @@ namespace Isis { 
      // Collect data to fit the polynomial over.

      std::vector<double> scaledTimes;

      std::vector< std::vector<double> > coordinates;

      // Compute how many points are required to fit the data

      //   This over estimates because the continuity conditions reduce the

      //   number of samples required.

      int requiredSamples = m_segments * ( p_degree + 1 ) * 3;

      // Double it to ensure that each segment is sufficiently supported

      requiredSamples *= 2;

      // If there are not enough samples, generate more

      if ( (int)p_cacheTime.size() < requiredSamples ) {

        double sampleRate = ( p_cacheTime.back() - p_cacheTime.front() ) / requiredSamples;

        for (int i = 0; i < requiredSamples; i++) {

          double sampleTime = p_cacheTime.front() + sampleRate * i;

          SetEphemerisTime( sampleTime );

          scaledTimes.push_back( (sampleTime - p_baseTime) / p_timeScale);

          coordinates.push_back( Coordinate() );

        }

      }

      else {

        for (int i = 0; i < (int)p_cacheTime.size(); i++) {

          SetEphemerisTime( p_cacheTime[i] );

          scaledTimes.push_back( (p_cacheTime[i] - p_baseTime) / p_timeScale);

          coordinates.push_back( Coordinate() );

        }

      }



      // Fit the polynomial.

      m_polynomial.fitPolynomials(scaledTimes, coordinates, m_segments);
@@ -988,6 +1006,81 @@ namespace Isis { 
  }





  /**

   * Create a test polynomial fit over the cached position data.

   * The polynomial works in scaled time, so all inputs must be adjusted

   * by subtracting the base time and then dividing by the time scale.

   * 

   * @return @b PiecewisePolynomial A piecewise polynomial fit over the cached

   *                                position data based on the segment count

   *                                and polynomial degree.

   */

  PiecewisePolynomial SpicePosition::testFit() {

    // Check to see if the position is already a Polynomial Function

    if (p_source == PolyFunction)

      return m_polynomial;



    // Check that there is data available to fit a polynomial to

    if ( p_cache.empty() ) {

      QString msg = "Cannot set a polynomial, no coordinate data is available.";

      throw IException(IException::Programmer, msg, _FILEINFO_);

    }



    // Adjust the degree of the polynomial to the available data

    int fitDegree = p_degree;

    if ( fitDegree > (int)p_cache.size() - 1 ) {

      fitDegree = p_cache.size() - 1;

    }



    // Recompute the time scaling

    ComputeBaseTime();



    // Compute first and last times in scaled time.

    //   If there is only one cached time, use the full extents.

    //   Otherwise scale the first and last times.

    double scaledFirstTime = -DBL_MAX;

    double scaledLastTime = DBL_MAX;

    if ( p_cacheTime.size() > 1) {

      scaledFirstTime = (p_cacheTime.front() - p_baseTime) / p_timeScale;

      scaledLastTime = (p_cacheTime.back() - p_baseTime) / p_timeScale;

    }



    // Initialize the zero polynomial.

    PiecewisePolynomial testPoly(scaledFirstTime, scaledLastTime, fitDegree, 3);



    // Collect data to fit the polynomial over.

    std::vector<double> scaledTimes;

    std::vector< std::vector<double> > coordinates;

    // Compute how many points are required

    //   This over estimates because the continuity conditions reduce the

    //   number of samples required.

    int requiredSamples = m_segments * ( p_degree + 1 ) * 3;

    // Double it to ensure that each segment is sufficiently supported

    requiredSamples *= 2;

    // If there are not enough samples, generate more

    if ( (int)p_cacheTime.size() < requiredSamples ) {

      double sampleRate = ( p_cacheTime.back() - p_cacheTime.front() ) / requiredSamples;

      for (int i = 0; i < requiredSamples; i++) {

        double sampleTime = p_cacheTime.front() + sampleRate * i;

        SetEphemerisTime( sampleTime );

        scaledTimes.push_back( (sampleTime - p_baseTime) / p_timeScale);

        coordinates.push_back( Coordinate() );

      }

    }

    else {

      for (int i = 0; i < (int)p_cacheTime.size(); i++) {

        SetEphemerisTime( p_cacheTime[i] );

        scaledTimes.push_back( (p_cacheTime[i] - p_baseTime) / p_timeScale);

        coordinates.push_back( Coordinate() );

      }

    }



    // Fit the polynomial.

    testPoly.fitPolynomials(scaledTimes, coordinates, m_segments);



    return testPoly;

  }





  /** 

   * Set the coefficients of a polynomial fit to

isis/src/base/objs/SpicePosition/SpicePosition.h

+2 −0
Original line number Diff line number Diff line
@@ -249,6 +249,8 @@ namespace Isis { 


      void SetPolynomial(const Source type = PolyFunction);



      PiecewisePolynomial testFit();



      void SetPolynomial(const std::vector<double>& XC,

                         const std::vector<double>& YC,

                         const std::vector<double>& ZC,