TgoCassisDistortionMap.cpp

Go to the documentation of this file.

#include <cmath>

#include <QDebug>

#include "IString.h"
#include "TgoCassisDistortionMap.h"

namespace Isis {
  TgoCassisDistortionMap::TgoCassisDistortionMap(Camera *parent, 
                                                 int naifIkCode) 
      : CameraDistortionMap(parent) {

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

    for(int i = 0; i < 6; i++) {
      m_A1_corr.push_back(p_camera->getDouble(od + "A1_CORR", i));
      m_A2_corr.push_back(p_camera->getDouble(od + "A2_CORR", i));
      m_A3_corr.push_back(p_camera->getDouble(od + "A3_CORR", i));
      m_A1_dist.push_back(p_camera->getDouble(od + "A1_DIST", i));
      m_A2_dist.push_back(p_camera->getDouble(od + "A2_DIST", i));
      m_A3_dist.push_back(p_camera->getDouble(od + "A3_DIST", i));
    }
    m_pixelPitch = p_camera->getDouble("INS" + toString(naifIkCode) + "_PIXEL_PITCH");
    m_width  = p_camera->getDouble("INS" + toString(naifIkCode) + "_FILTER_SAMPLES");
    m_height = p_camera->getDouble("INS" + toString(naifIkCode) + "_FILTER_LINES");
  }


  TgoCassisDistortionMap::~TgoCassisDistortionMap() {
  }


  bool TgoCassisDistortionMap::SetFocalPlane(const double dx, 
                                             const double dy) {

    p_focalPlaneX = dx;
    p_focalPlaneY = dy;

    // calculate divisor using A3_corr coeffiecients
    double divider = chiDotA(dx, dy, m_A3_corr);
    // This distortion model is only valid for values on the CCD:
    // -1/2 * pixel pitch * CCD width  = -10.24 < x < 10.24 = 1/2 * pixel pitch * CCD width
    // -1/2 * pixel pitch * CCD height = -10.24 < y < 10.24 = 1/2 * pixel pitch * CCD height
    // 
    // Also, the zeros for the divider variable fall well outside the boundary
    // of the CCD. (See $ISIS3DATA/tgo/assets/distortion/DistortionModelA3CorrRoots.jpg).
    // 
    // So, whenever x or y are too far from center or divider is near zero,
    // return the given inputs
    
    if ( dx < -0.5*m_pixelPitch*m_width  - 0.2 ||
         dx >  0.5*m_pixelPitch*m_width  + 0.2 ||
         dy < -0.5*m_pixelPitch*m_height - 0.2 ||
         dy >  0.5*m_pixelPitch*m_height + 0.2 ) {
      p_undistortedFocalPlaneX = dx;
      p_undistortedFocalPlaneY = dy;

      return true;
    }

    // get undistorted ideal (x,y) coordinates
    double ux = chiDotA(dx, dy, m_A1_corr) / divider;
    double uy = chiDotA(dx, dy, m_A2_corr) / divider;

    p_undistortedFocalPlaneX = ux;
    p_undistortedFocalPlaneY = uy;

    return true;
  }


  bool TgoCassisDistortionMap::SetUndistortedFocalPlane(const double ux,
                                                        const double uy) {
    p_undistortedFocalPlaneX = ux;
    p_undistortedFocalPlaneY = uy;

    // calculate divisor using A3_dist coeffiecients
    double divider = chiDotA(ux, uy, m_A3_dist);
    // This distortion model is only valid for values on the CCD:
    // -1/2 * pixel pitch * CCD width  = -10.24 < x < 10.24 = 1/2 * pixel pitch * CCD width
    // -1/2 * pixel pitch * CCD height = -10.24 < y < 10.24 = 1/2 * pixel pitch * CCD height
    // 
    // Also, the zeros for the divider variable fall well outside the boundary
    // of the CCD. (See $ISIS3DATA/tgo/assets/distortion/DistortionModelA3DistRoots.jpg).
    // 
    // So, whenever x or y are too far from center or divider is near zero,
    // return the given inputs

    if ( ux < -0.5*m_pixelPitch*m_width  - 0.2 ||
         ux >  0.5*m_pixelPitch*m_width  + 0.2 ||
         uy < -0.5*m_pixelPitch*m_height - 0.2 ||
         uy >  0.5*m_pixelPitch*m_height + 0.2 ) {
      p_focalPlaneX = ux;
      p_focalPlaneY = uy;

      return true;
    }

    // get undistorted ideal (x,y) coordinates
    double dx = chiDotA(ux, uy, m_A1_dist) / divider;
    double dy = chiDotA(ux, uy, m_A2_dist) / divider;

    p_focalPlaneX = dx;
    p_focalPlaneY = dy;
    return true;
  }


  double TgoCassisDistortionMap::chiDotA(double x,
                                         double y,
                                         QList<double> A) {  
    double x2 = x*x;
    double y2 = y*y;
    double xy = x*y;

    return A[0] * x2
         + A[1] * xy 
         + A[2] * y2 
         + A[3] * x 
         + A[4] * y 
         + A[5];
  }
}