TgoCassisDistortionMap.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#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 $ISISDATA/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 $ISISDATA/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];
  }
}