Interpolator.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include <string>
#include "IException.h"
#include "Interpolator.h"
 
using namespace std;
namespace Isis {
  Interpolator::Interpolator() {
    Init();
  }

  Interpolator::Interpolator(const interpType &type) {
    Init();
    p_type = type;
  }

  Interpolator::~Interpolator() {}

  void Interpolator::Init() {
    p_type = None;
  }

  double Interpolator::Interpolate(const double isamp, const double iline,
                                   const double buf[]) {
 
    switch(p_type) {
      case None: {
          string message = "Interpolator type not set";
          throw IException(IException::Programmer, message, _FILEINFO_);
          break;
        }
      case NearestNeighborType:
        return NearestNeighbor(isamp, iline, buf);
        break;
      case BiLinearType:
        return BiLinear(isamp, iline, buf);
        break;
      case CubicConvolutionType:
        return CubicConvolution(isamp, iline, buf);
        break;
    }
 
    string message = "Invalid interpolator";
    throw IException(IException::Programmer, message, _FILEINFO_);
  }

  void Interpolator::SetType(const interpType &type) {
    p_type = type;
  }
 

  double Interpolator::NearestNeighbor(const double isamp, const double iline,
                                       const double buf[]) {
    return buf[0];
  }
 

  double Interpolator::BiLinear(const double isamp, const double iline,
                                const double buf[]) {
 
    // Get the fractional portions of the sample and line coordinates
    double j = int (isamp);
    double k = int (iline);
    double a = isamp - j;
    double b = iline - k;
 
    // If any of the 4 pixels are special pixels, drop down to a nearest neighbor
    for(int i = 0; i < 4; i++) {
      if(Isis::IsSpecial(buf[i])) {
        return NearestNeighbor(isamp, iline,
                               &buf[(int)(a + 0.5) + 2 * (int)(b+0.5)]);
        break;
      }
    }
 
    // Otherwise do the bilinear
    return (1.0 - a) * (1.0 - b) * buf[0] +
           a * (1.0 - b) * buf[1] +
           (1.0 - a) * b * buf[2] +
           a * b * buf[3];
 
  }

  double Interpolator::CubicConvolution(const double isamp,
                                        const double iline,
                                        const double buf[]) {
 
    // If any of the 16 pixels are special pixels, drop down to a bilinear
    for(int i = 0; i < 16; i++) {
      if(Isis::IsSpecial(buf[i])) {
        double tbuf[4] = {buf[5], buf[6], buf[9], buf[10]};
        return BiLinear(isamp, iline, tbuf);
        break;
      }
    }
 
    double j = int (isamp);
    double k = int (iline);
    double a = isamp - j;
    double b = iline - k;

    return  -b * (1.0 - b) * (1.0 - b) * (-a * (1.0 - a) * (1.0 - a) * buf[0] +
                                          (1.0 - 2.0 * a * a + a * a * a) * buf[1] +
                                          a * (1.0 + a - a * a) * buf[2] -
                                          a * a * (1.0 - a) * buf[3]) +
 
            (1.0 - 2.0 * b * b + b * b * b) * (-a * (1.0 - a) * (1.0 - a) * buf[4] +
                (1.0 - 2.0 * a * a + a * a * a) * buf[5] +
                a * (1.0 + a - a * a) * buf[6] -
                a * a * (1.0 - a) * buf[7]) +
 
            b * (1.0 + b - b * b) * (-a * (1.0 - a) * (1.0 - a) * buf[8] +
                                     (1.0 - 2.0 * a * a + a * a * a) * buf[9] +
                                     a * (1.0 + a - a * a) * buf[10] -
                                     a * a * (1.0 - a) * buf[11]) +
 
            b * b * (b - 1.0) * (-a * (1.0 - a) * (1.0 - a) * buf[12] +
                                 (1.0 - 2.0 * a * a + a * a * a) * buf[13] +
                                 a * (1.0 + a - a * a) * buf[14] -
                                 a * a * (1.0 - a) * buf[15]);
 
  }

  int Interpolator::Samples() {
 
    switch(p_type) {
      case None: {
          string message = "Interpolator type not set";
          throw IException(IException::Programmer, message, _FILEINFO_);
          break;
        }
      case NearestNeighborType:
        return 1;
        break;
      case BiLinearType:
        return 2;
        break;
      case CubicConvolutionType:
        return 4;
        break;
    }
 
    string message = "Invalid interpolator";
    throw IException(IException::Programmer, message, _FILEINFO_);
  }

  int Interpolator::Lines() {
 
    switch(p_type) {
      case None: {
          string message = "Interpolator type not set";
          throw IException(IException::Programmer, message, _FILEINFO_);
          break;
        }
      case NearestNeighborType:
        return 1;
        break;
      case BiLinearType:
        return 2;
        break;
      case CubicConvolutionType:
        return 4;
        break;
    }
 
    string message = "Invalid interpolator";
    throw IException(IException::Programmer, message, _FILEINFO_);
  }

  double Interpolator::HotSample() {
 
    switch(p_type) {
      case None: {
          string message = "Interpolator type not set";
          throw IException(IException::Programmer, message, _FILEINFO_);
          break;
        }
        // To get the correct pixel for NN you have to round the sample
      case NearestNeighborType:
        return -0.5;
        break;
        // To get the correct pixel for BL you have to truncate the sample
      case BiLinearType:
        return 0.0;
        break;
        // To get the correct pixel for CC you have to truncate the sample
      case CubicConvolutionType:
        return 1.0;
        break;
    }
 
    string message = "Invalid interpolator";
    throw IException(IException::Programmer, message, _FILEINFO_);
  }

  double Interpolator::HotLine() {
 
    switch(p_type) {
      case None: {
          string message = "Interpolator type not set";
          throw IException(IException::Programmer, message, _FILEINFO_);
          break;
        }
        // To get the correct pixel for NN you have to round the line
      case NearestNeighborType:
        return -0.5;
        break;
        // To get the correct pixel for BL you have to truncate the line
      case BiLinearType:
        return 0.0;
        break;
        // To get the correct pixel for CC you have to truncate the line
      case CubicConvolutionType:
        return 1.0;
        break;
    }
 
    string message = "Invalid interpolator";
    throw IException(IException::Programmer, message, _FILEINFO_);
  }
}