Commit bca4f132 authored by Jeannie Backer's avatar Jeannie Backer
Browse files

Added distortion model to TGO CaSSIS camera model. References #4593 

git-svn-id: http://subversion.wr.usgs.gov/repos/prog/isis3/trunk@7520 41f8697f-d340-4b68-9986-7bafba869bb8
parent 217e3a87

isis/src/tgo/objs/TgoCassisCamera/TgoCassisCamera.cpp

+3 −3
Original line number Diff line number Diff line
@@ -19,6 +19,7 @@ 
 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on

 *   http://www.usgs.gov/privacy.html.

 */

#include "TgoCassisCamera.h"



#include <cmath>
@@ -26,9 +27,7 @@ 
#include <QString>

#include <QVariant>



#include "TgoCassisCamera.h"

#include "CameraDetectorMap.h"

#include "CameraDistortionMap.h"

#include "CameraFocalPlaneMap.h"

#include "CameraGroundMap.h"

#include "CameraSkyMap.h"
@@ -37,6 +36,7 @@ 
#include "IString.h"

#include "iTime.h"

#include "NaifStatus.h"

#include "TgoCassisDistortionMap.h"



using namespace std;
@@ -111,7 +111,7 @@ namespace Isis { 
    detMap->SetStartingDetectorLine(getDouble("INS" + filter + "_FILTER_OFFSET"));



    // Setup distortion map

    new CameraDistortionMap(this);

    new TgoCassisDistortionMap(this, naifIkCode());



    // Setup the ground and sky map

    new CameraGroundMap(this);

isis/src/tgo/objs/TgoCassisCamera/TgoCassisCamera.h

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

   * </pre>

   *  

   * @ingroup SpiceInstrumentsAndCameras

   * @ingroup Tgo

   * @author 2017-01-26 Kris Becker

   *

   * @internal

isis/src/tgo/objs/TgoCassisCamera/TgoCassisCamera.truth

+17 −17
Original line number Diff line number Diff line 
Unit Test for TgoCassisCamera...

FileName: CAS-MCO-2016-11-22T16.38.39.354-NIR-02036-B1.cub

FileName:  "CAS-MCO-2016-11-22T16.38.39.354-NIR-02036-B1.cub"

CK Frame:  -143420



Kernel IDs:
@@ -10,23 +10,23 @@ SPK Reference ID = 1 


Shutter open =  533104787.536896408

Shutter close =  533104787.538624406

Focal Length = 880.000000000

Focal Length =  880



For upper left corner ...

DeltaSample = 0.000000000

DeltaLine = 0.000000000

DeltaSample =  0

DeltaLine =  0



For upper right corner ...

DeltaSample = 0.000000000

DeltaLine = 0.000000000

DeltaSample =  0

DeltaLine =  0



For lower left corner ...

DeltaSample = 0.000000000

DeltaLine = 0.000000000

DeltaSample =  0

DeltaLine =  0



For lower right corner ...

DeltaSample = 0.000000000

DeltaLine = 0.000000000

DeltaSample =  0

DeltaLine =  0



For center pixel position ...

Latitude OK

isis/src/tgo/objs/TgoCassisCamera/TgoCassisDistortionMap.cpp

0 → 100644
+249 −0
Original line number Diff line number Diff line 
/**

 * @file

 * $Revision: 1.2 $

 * $Date: 2008/11/24 16:40:30 $

 *

 *   Unless noted otherwise, the portions of Isis written by the USGS are

 *   public domain. See individual third-party library and package descriptions

 *   for intellectual property information, user agreements, and related

 *   information.

 *

 *   Although Isis has been used by the USGS, no warranty, expressed or

 *   implied, is made by the USGS as to the accuracy and functioning of such

 *   software and related material nor shall the fact of distribution

 *   constitute any such warranty, and no responsibility is assumed by the

 *   USGS in connection therewith.

 *

 *   For additional information, launch

 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html

 *   in a browser or see the Privacy &amp; Disclaimers page on the Isis website,

 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on

 *   http://www.usgs.gov/privacy.html.

 */

#include <cmath>



#include <QDebug>



#include "IString.h"

#include "TgoCassisDistortionMap.h"



namespace Isis {

  /** Camera distortion map constructor

   *

   * Create a camera distortion map.  This class maps between distorted

   * and undistorted focal plane x/y's.  The default mapping is the

   * identity, that is, the focal plane x/y and undistorted focal plane

   * x/y will be identical.

   *

   * @param parent        the parent camera that will use this distortion map

   * @param zDirection    the direction of the focal plane Z-axis

   *                      (either 1 or -1)

   *

   */

  TgoCassisDistortionMap::TgoCassisDistortionMap(Camera *parent, 

                                                 int naifIkCode) 

      : CameraDistortionMap(parent) {



    QString od = "INS" + toString(naifIkCode) + "_OD_";



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

      m_A1.push_back(p_camera->getDouble(od + "A1", i));

      m_A2.push_back(p_camera->getDouble(od + "A2", i));

      m_A3.push_back(p_camera->getDouble(od + "A3", i));

    }

  }





  /** Compute undistorted focal plane x/y

   *

   * Compute undistorted focal plane x/y given a distorted focal plane x/y.

   *

   * @param dx distorted focal plane x in millimeters

   * @param dy distorted focal plane y in millimeters

   *

   * @return if the conversion was successful

   * @see SetDistortion

   * @todo Generalize polynomial equation

   */

  bool TgoCassisDistortionMap::SetFocalPlane(const double dx, 

                                             const double dy) {



    p_focalPlaneX = dx;

    p_focalPlaneY = dy;



    // i and j are normalized distorted coordinates

    double i = normalize(dx);

    double j = normalize(dy);



    // convenience variables

    double i2 = i*i;

    double j2 = j*j;

    double ij = i*j;



    // divider we might nead for debuging

    double divider = i2 * m_A3[0] 

                     + ij * m_A3[1] 

                     + j2 * m_A3[2] 

                     + i  * m_A3[3] 

                     + j  * m_A3[4] 

                     + 1;



    double xNorm = ( i2 * m_A1[0] 

                     + ij * m_A1[1] 

                     + j2 * m_A1[2] 

                     + i  * m_A1[3] 

                     + j  * m_A1[4] 

                     + m_A1[5] )

                   / divider;



    double yNorm = ( i2 * m_A2[0] 

                     + ij * m_A2[1] 

                     + j2 * m_A2[2] 

                     + i  * m_A2[3] 

                     + j  * m_A2[4] 

                     + m_A2[5] )

                   / divider;



    // denormalize ideal (x,y) coordinates

    p_undistortedFocalPlaneX = denormalize(xNorm);

    p_undistortedFocalPlaneY = denormalize(yNorm);



    return true;

  }





  /** Compute distorted focal plane x/y

   *

   * Compute distorted focal plane x/y given an undistorted focal plane x/y.

   *

   * @param ux undistorted focal plane x in millimeters

   * @param uy undistorted focal plane y in millimeters

   *

   * @return if the conversion was successful

   * @see SetDistortion

   * @todo Generalize polynomial equation

   * @todo Figure out a better solution for divergence condition

   */

  bool TgoCassisDistortionMap::SetUndistortedFocalPlane(const double ux,

                                                        const double uy) {

    p_undistortedFocalPlaneX = ux;

    p_undistortedFocalPlaneY = uy;



    // x, y are normalized undistorted (ideal) coordinates

    double xNorm = normalize(ux);

    double yNorm = normalize(uy);



    // i, j are distorted coordinates

    double iNorm = xNorm;

    double jNorm = yNorm;

    int newtonIterations = 2000;  //???

    // newton's method to iterate

    for( int index = 0; index < newtonIterations; ++index ) {

      /* 

       * compute F(i_pre, j_pre)

       * 

       * F_vec(i,j) = [ initialX - ( A1 * chi' ) / ( A3 * chi' ) ]

       *              [ initialY - ( A2 * chi' ) / ( A3 * chi' ) ]

      */



      if (SetFocalPlane(iNorm, jNorm) ) {

        double xPredict = p_undistortedFocalPlaneX;

        double yPredict = p_undistortedFocalPlaneY;

        double divider = iNorm * iNorm * m_A3[0] 

                         + iNorm * jNorm * m_A3[1] 

                         + jNorm * jNorm * m_A3[2] 

                         + iNorm * m_A3[3] 

                         + jNorm * m_A3[4] 

                         + 1;



        double f11 = xNorm - xPredict;

        double f21 = yNorm - yPredict;



        /* 

         * compute the Jacobian, J(i, j):

         * 

         * J(i,j) = [ - ( (A1 * dchi / di') - xPredicted*(A3 * dchi / di) ) / (A3 * chi'),...

         *            - ( (A1 * dchi / dj') - xPredicted*(A3 * dchi / dj) ) / (A3 * chi');...

         *            - ( (A2 * dchi / di') - yPredicted*(A3 * dchi / di) ) / (A3 * chi'),...

         *            - ( (A2 * dchi / dj') - yPredicted*(A3 * dchi / dj) ) / (A3 * chi')]

         * 

         * dchi / di = [2i   j   0   1   0   0]

         * dchi / dj = [ 0   i  2j   0   1   0]

         */

        double j11 = - ( ( 2 * iNorm * m_A1[0] + jNorm * m_A1[1] + m_A1[3] ) 

                         - xPredict * ( 2 * iNorm * m_A3[0] + jNorm * m_A3[1] + m_A3[3] ) ) 

                       / divider;

        double j12 = - ( ( iNorm * m_A1[1] + 2 * jNorm * m_A1[2] + m_A1[4] ) 

                         - xPredict * ( iNorm * m_A3[1] + 2 * jNorm * m_A3[2] + m_A3[4] ) )

                       / divider;

        double j21 = - ( ( 2 * iNorm * m_A2[0] + jNorm * m_A2[1] + m_A2[3] ) 

                         - yPredict * ( 2 * iNorm * m_A3[0] + jNorm * m_A3[1] + m_A3[3] ) )

                       / divider;

        double j22 = - ( ( iNorm * m_A2[1] + 2 * jNorm * m_A2[2] + m_A2[4] ) 

                         - yPredict * ( iNorm * m_A3[1] + 2 * jNorm * m_A3[2] + m_A3[4] ) )

                       / divider;



        /* 

         * compute update

         * 

         * [i_n, j_n]  = [i_n-1, j_n-1] - inv(J(i_n-1, j_n-1))*F(i_n-1, j_n-1);

         * 

         *                    1             [  J22 -J12 ] [ F11 ]   [  J22*F11 - J12*F21 ]

         * inv(J)*F =  ------------------ * [ -J21  J11 ] [ F21 ] = [ -J21*F11 + J11*F21 ]

         *             J11*J22 - J12*J21

         */

        double di = - ( j22*f11 - j12*f21) / (j11*j22 - j12*j21);

        double dj = - (-j21*f11 + j11*f21) / (j11*j22 - j12*j21);

        iNorm = iNorm + di;

        jNorm = jNorm + dj;

      }

      else {

        return false;

      }



    }

    p_undistortedFocalPlaneX = ux;

    p_undistortedFocalPlaneY = uy;

    // denormalize distorted (i,j) coordinates

    p_focalPlaneX = denormalize(iNorm);

    p_focalPlaneY = denormalize(jNorm);

    return true;

  }





  /**

   * Normalize the value using the dimensions of the CCD, 2048 x 2048. This 

   * method uses the formula below to normalize the given value: 

   *  

   *  @f[ norm = \frac{ value - \frac{2048}{2} }{ 2048 + 2048 } @f]

   * 

   * @param value Value to be normalize.

   * 

   * @return @b double The normalized value.

   */

  double TgoCassisDistortionMap::normalize(double value) {

    // use scaling factor based on size of CCD:

    // CCD width = CCD height = 2048

    return (value - 2048 / 2) / 4096;

  }





  /**

   * De-normalize the value using the dimensions of the CCD, 2048 x 2048. 

   * This method is the inverse function of the normalize() function. It 

   * uses the formula below to denormalize the given value: 

   *  

   *  @f[ denorm = value * (2048 + 2048) + \frac{2048}{2} @f]

   * 

   * @param value Value to be normalize.

   * 

   * @return @b double The normalized value.

   */

  double TgoCassisDistortionMap::denormalize(double value, double a, double b) {

    // use scaling factor based on size of CCD:

    // CCD width = CCD height = 2048

    return value * 4096 + 2048 / 2;

  }

  

}

isis/src/tgo/objs/TgoCassisCamera/TgoCassisDistortionMap.h

0 → 100644
+67 −0
Original line number Diff line number Diff line 
#ifndef TgoCassisDistortionMap_h

#define TgoCassisDistortionMap_h

/**

 * @file

 * $Revision: 1.2 $

 * $Date: 2008/11/24 16:40:31 $

 *

 *   Unless noted otherwise, the portions of Isis written by the USGS are

 *   public domain. See individual third-party library and package descriptions

 *   for intellectual property information, user agreements, and related

 *   information.

 *

 *   Although Isis has been used by the USGS, no warranty, expressed or

 *   implied, is made by the USGS as to the accuracy and functioning of such

 *   software and related material nor shall the fact of distribution

 *   constitute any such warranty, and no responsibility is assumed by the

 *   USGS in connection therewith.

 *

 *   For additional information, launch

 *   $ISISROOT/doc//documents/Disclaimers/Disclaimers.html

 *   in a browser or see the Privacy &amp; Disclaimers page on the Isis website,

 *   http://isis.astrogeology.usgs.gov, and the USGS privacy and disclaimers on

 *   http://www.usgs.gov/privacy.html.

 */



#include <vector>

#include "CameraDistortionMap.h"



namespace Isis {

  /** 

   *  Distort/undistort focal plane coordinates

   *

   * Creates a map for adding/removing optical distortions

   * from the focal plane of a camera.

   *

   * @ingroup SpiceInstrumentsAndCameras

   * @ingroup Tgo

   *

   * @see TgoCassisCamera

   *

   * @author 2008-08-22 Steven Lambright

   *

   * @internal

   *   @history 2017-04-03 Jeannie Walldren - Original version.

   */

  class TgoCassisDistortionMap : public CameraDistortionMap {

    public:

      TgoCassisDistortionMap(Camera *parent, int naifIkCode);



      //! Destructor

      virtual ~TgoCassisDistortionMap() {};



      virtual bool SetFocalPlane(const double dx, const double dy);



      virtual bool SetUndistortedFocalPlane(const double ux, const double uy);



    private:

      double normalize(double value, double a, double b);

      double denormalize(double value, double a, double b);

      double divisor(double i, double j);



      QList<double> m_A1; //!< First row of parameters of rational distortion model.

      QList<double> m_A2; //!< Second row of parameters of rational distortion model.

      QList<double> m_A3; //!< Third row of parameters of rational distortion model.

  };

};

#endif