Cleaned up hyb2onccal to match the JAXA paper (#3952)

      process FIT files using the old keywords.  Also set TARGET=RYUGU in the event the user

      does not specify a target, and the target in the fit label is set to SKY.

    </change>

    <change name="Kaitlyn Lee" date="2020-07-21">

      Updated calibration method to match JAXAs calibration by:

      1. moving linearity correction to after smear correction,

      2. commenting out dark current and linearity correction,

      3. changing an addition to subtraction in the g_bias calculation,

      4. refactoring the smear correction to be before the main loop to speed up the program

         and to be dependent on the bit-corrected input image DNs instead of the already corrected

         output image DNs,

      5. and removing `*g_sensitivity*g_texp` in the flat-field correction so that the output

         aligns with l2b data.

    </change>

  </history>

  <category>

#ifndef Hyb2OncCalUtils_h

#define Hyb2OncCalUtils_h

#include "CSVReader.h"

#include "IException.h"

#include "LineManager.h"

#include "NaifStatus.h"

#include <string>

#include <vector>

#include <algorithm>

#include "AlphaCube.h"

#include "Buffer.h"

#include "CSVReader.h"

#include "FileName.h"

#include "IException.h"

#include "iTime.h"

#include "LineManager.h"

#include "NaifStatus.h"

#include "Pixel.h"

#include "ProcessByLine.h"

#include "ProcessBySample.h"

#include "Statistics.h"

#include "TextFile.h"

// OpenCV libraries

#include <opencv2/opencv.hpp>

/**

 * @author 2016-04-04 Tyler Wilson

 *

 *

 */

using namespace cv;

using namespace std;

static QString g_iofCorrection("IOF");  //!< Is I/F correction to be applied?

// I/F variables

//static double g_iof(1.0);

static double g_iofScale(1.0);

static double g_solarFlux(1.0);  //!< The solar flux (used to calculate g_iof).

static double g_sensitivity(1.0);

static double g_effectiveBandwidth(1.0);

static double g_J(1.0);

namespace Isis {

// Temporary cube file pointer deleter

struct TemporaryCubeDeleter {

  static inline void cleanup(Cube *cube) {

    if (cube) {

      FileName filename(cube->fileName());

      delete cube;

      remove(filename.expanded().toLatin1().data());

InstrumentType g_instrument;

static AlphaCube *alpha(0);

static Pvl g_configFile;

 * @return The value of the function at the point x.

 */

double linearFun(double Iobs,double x, double g[3]) {

  return Iobs - g[0]*x -g[1]*pow(x,2.0) -g[2]*pow(x,3.0);

  return Iobs - (g[0] * x) - (g[1] * pow(x, 2.0)) - (g[2] * pow(x, 3.0));

}

/**

 * @return

 */

double dFun(double x, double g[3]) {

  return -g[0] - 2*g[1]*x -3*g[2]*pow(x,2.0);

  return -g[0] - (2 * g[1] * x) - (3 * g[2] * pow(x, 2.0));

}

   int iter = 0;

   int maxIterations = 500;

   x[0] = x0;

   while (dx > epsilon)  {

   while (dx > epsilon) {

     x[1] = x[0] - linearFun(Iobs, x[0], g) / dFun(x[0], g);

     dx = fabs(x[1] - x[0]);

     x[0] = x[1];

     iter++;

     if (iter > maxIterations) {

       return false;

     }

   }

*   @history 2019-02-12 Tyler Wilson - Modified to support new calibration settings/formulas.

*/

void Calibrate(vector<Buffer *>& in, vector<Buffer *>& out) {

  Buffer& imageIn   = *in[0];

  Buffer& flatField = *in[1];

  Buffer& imageOut  = *out[0];

  int alphaSample = alpha->AlphaSample(currentSample);

  if ( (alphaSample <= pixelsToNull) || (alphaSample >= (1024 - pixelsToNull)) ) {

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

      imageOut[i] = Isis::Null;

    }

    return;

  }

  double smear = 0;

  for (int j = 0; j < imageIn.size(); j++) {

    // Left out g_darkCurrent subtraction for now.

    smear += ((imageIn[j] * pow(2.0, 12 - g_bitDepth) - g_bias) / imageIn.size());

  }

  smear *= g_timeRatio;

  // Iterate over the line space

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

    imageOut[i] = imageIn[i] * pow(2.0, 12 - g_bitDepth);

    // Apply compression factor here to raise LOSSY dns to proper response

    // 1) BIAS Removal - Only needed if not on-board corrected

    // BIAS Removal - Only needed if not on-board corrected

    if (!g_onBoardSmearCorrection) {

      if ( (imageOut[i] - g_bias) <= 0.0) {

        imageOut[i] = Null;

        continue;

      }

    }

    double dn = imageOut[i];    

    double result = 1.0;

    double x0 = 1.0;

    newton_rapheson(dn,x0, g_L,result );   

    imageOut[i] = result;

    // DARK Current

    imageOut[i] = imageOut[i] - g_darkCurrent;    

    // Smear correction

    if (!g_onBoardSmearCorrection) {

      double smear = 0;

      for (int j=0;j < imageIn.size();j++) {

        smear += (imageOut[j]/imageIn.size() );

      }

      smear*=g_timeRatio;

      imageOut[i] = imageOut[i] - smear;

    }

    // Both dark current and linearity correction are commented out for now since the JAXA team

    // does not currently do these steps.

    // DARK Current

    // imageOut[i] = imageOut[i] - g_darkCurrent;    

    // Linearity Correction

    // double dn = imageOut[i];    

    // double result = 1.0;

    // double x0 = 1.0;

    // newton_rapheson(dn,x0, g_L,result );   

    // imageOut[i] = result;

    // FLATFIELD correction

    //  Check for any special pixels in the flat field (unlikely)

      }

      else {

        if (flatField[i] != 0) {

          imageOut[i] /= (flatField[i]*g_sensitivity*g_texp);

          imageOut[i] /= (flatField[i]);

        }

      }

    }

  }

  return;

}

/**

 * @brief Translates a 1-banded Isis::Cube to an OpenMat object

 *

 * @return @b Mat A pointer to the OpenMat object

 */

Mat *isis2mat(Cube *icube) {

  int nlines = icube->lineCount();

  int nsamples = icube->sampleCount();

  Mat *matrix = new Mat(nlines, nsamples, CV_64F);

  // Set up line manager and read in the data

  LineManager linereader(*icube);

  for (int line = 0; line < nlines; line++) {

      matrix->at<double>(line, samp) = (double)linereader[samp];

    }

 }

return matrix;

}

 * @param matrix A pointer to the OpenMat object

 *

 * @param cubeName The name of the Isis::Cube that is being created.

 *

 */

void mat2isis(Mat *matrix, QString cubeName) {

    for ( int samp=0; samp < nsamples; samp++ ) {

      linewriter[samp] = matrix->at<double>(linewriter.Line() - 1,samp);

    }

    ocube.write(linewriter);

  }

  ocube.close();

}

 * @param matrix A pointer to the OpenMat object

 *

 * @param cubeName The name of the ISIS::Cube that is being created.

 *

 */

void translate(Cube *flatField,double *transform, QString fname) {

  Mat *originalMat = isis2mat(flatField);

  double scale = transform[0];

  Mat temp = *originalMat;

  Mat originalCropped = temp(Rect(startsample, startline, width, height));

  if (scale == 1) {

*   @history 2019-02-12 Tyler Wilson - Modified to support new calibration settings/formulas.

*/

FileName DetermineFlatFieldFile(const QString &filter) {

  QString fileName = "$hayabusa2/calibration/flatfield/";

  // FileName consists of binned/notbinned, camera, and filter

  fileName += "flat_" + filter.toLower() + "_norm.cub";

  FileName final(fileName);

  return final;

}

* @param config QString Name of the calibration file to load.

*

* Loads g_b0-g_b2,g_bae0-g_bae1,g_d0-g_d1

*

*

*/

QString loadCalibrationVariables(const QString &config)  {

  //  UserInterface& ui = Application::GetUserInterface();

  FileName calibFile(config);

  if ( config.contains("?") ) calibFile = calibFile.highestVersion();

  // Load Smear Removal Variables

  g_Tvct = Smear["Tvct"];

  // Load DarkCurrent variables and calculate the dark current

  g_d0 = DarkCurrent["D"][0].toDouble();

  g_d1 = DarkCurrent["D"][1].toDouble();

      CCDTemp = g_CCD_T_temperature;

  }

   g_darkCurrent = g_texp*exp(0.10*CCDTemp+0.52);

   g_darkCurrent = g_texp * exp(g_d0 * (CCDTemp + g_d1));

  // Load Bias variables

  g_b0 = Bias["B"][0].toDouble();

  g_bae0 = Bias["B_AE"][0].toDouble();

  g_bae1 = Bias["B_AE"][1].toDouble();

  // Compute BIAS correction factor (it's a constant so do it once!)

  g_bias = g_b0+g_b1*g_CCD_T_temperature+g_b2*g_ECT_T_temperature;

  g_bias = g_b0 + (g_b1 * g_CCD_T_temperature) + (g_b2 * g_ECT_T_temperature);

  g_bias *= (g_bae0 + g_bae1*g_AEtemperature); //bias correction factor

  g_bias *= (g_bae0 - (g_bae1 * g_AEtemperature)); //bias correction factor

  // Load the Solar Flux for the specific filter

  g_solarFlux = solar[g_filter.toLower()];

  g_J = g_solarFlux / (g_effectiveBandwidth * .0001);

  // Load the linearity variables

  g_L[0] = linearity["L"][0].toDouble();

  g_L[1] = linearity["L"][1].toDouble();

  g_L[2] = linearity["L"][2].toDouble();

  return ( calibFile.original() );

  return calibFile.original();

}

#include "Statistics.h"

#include "TextFile.h"

using namespace std;

namespace Isis {

  void hyb2onccal(UserInterface &ui) {

    // g_iofCorrection = ui.GetString("UNITS");

    const QString hyb2cal_program = "hyb2onccal";

    QString hyb2cal_runtime = Application::DateTime();

    ProcessBySample p;

    Cube *icube = p.SetInputCube("FROM");

    // Basic assurances...

    try {

      g_filter = bandbin["FilterName"][0];

    }

    catch(IException &e) {

      QString msg = "Unable to read FilterName keyword in the BandBin group "

                    "from input file [" + icube->fileName() + "]";

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

    }

    QString instrument("");

      QString msg = "Unable to read InstrumentId keyword in the Instrument group "

                    "from input file [" + icube->fileName() + "]";

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

    }

    if (instrument=="ONC-W1") {

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

    }

    //Set up binning and image subarea mapping

    binning = inst["Binning"];

    int startLine = inst["SelectedImageAreaY1"];

    }

    if (g_bitDepth < 0) {

      g_bitDepth = 12;  //Correpsonds to no correction being done for bit depth

    }

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

    }

    try {

      g_AEtemperature = inst["ONCAETemperature"][0].toDouble();

    }

    catch(IException &e) {

      QString msg = "Unable to read [ONCAETemperature] keyword in the Instrument group "

                    "from input file [" + icube->fileName() + "]";

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

    }

    try {

      g_CCD_T_temperature = inst["ONCTCCDTemperature"][0].toDouble();

    }

    catch(IException &e) {

      QString msg = "Unable to read [ONCTCCDTemperature] keyword in the Instrument group "

    try {

      g_ECT_T_temperature = inst["ONCTElectricCircuitTemperature"][0].toDouble();

    }

    catch(IException &e) {

      QString msg = "Unable to read [ONCTElectricCircuitTemperature] keyword in the Instrument group "

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

    }

    try {

      g_startTime = inst["SpacecraftClockStartCount"][0];

    }

    catch (IException &e) {

      QString msg = "Unable to read [SpacecraftClockStartCount] keyword in the Instrument group "

                    "from input file [" + icube->fileName() + "]";

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

    }

    QString smearCorrection;

    try {

      smearCorrection = inst["SmearCorrection"][0];

    }

    catch (IException &e) {

      QString msg = "Unable to read [SmearCorrection] keyword in the Instrument group "

                    "from input file [" + icube->fileName() + "]";

    }

    else {

      qDebug() << icube->fileName();

    }

    QString compmode = inst["Compression"];

    // TODO: verify that the compression factor/scale is actually 16 for compressed Hayabusa2 images.

    g_compfactor = ("lossy" == compmode.toLower()) ? 16.0 : 1.0;

    QString target = inst["TargetName"];

    g_target = target;

    // NOTE we do not have a valid flat-field for the W1 or W2 images.

    FileName flatfile = "NONE";

    PvlGroup alphaCube;

      // Image is not cropped

      if (!g_cropped) {

        // Determine if we need to subsample the flat field if pixel binning occurred

        // TODO: test a binned image (add test case).

        if (binning > 1) {

            remove(reducedFlat.toLatin1().data());

            throw ie;

          }

        //QScopedPointer<Cube, TemporaryCubeDeleter> reduced(new Cube(reducedFlat, "r"));

        //flatcube.swap( reduced );

        }

        // Set up processing for flat field as a second input file

        CubeAttributeInput att;

        p.SetInputCube(reducedFlat, att);

      }

      // The image was cropped, so pull the same subarea from the flat file into a temp file

      else {

      // Image is cropped so we have to deal with it

        FileName transFlat =

        FileName::createTempFile("$TEMPORARY/" + flatfile.baseName() + "_translated.cub");

        Cube *flatOriginal = new Cube(flatfile.expanded() );

        double alphaStartSample = alphaCube["AlphaStartingSample"][0].toDouble();

        double alphaStartLine = alphaCube["AlphaStartingLine"][0].toDouble();

        double alphaEndSample = alphaCube["AlphaEndingSample"][0].toDouble();

        // Translates and scales the flatfield image.  Scaling

        // might be necessary in the event that the raw image was also binned.

        translate(flatOriginal,transform,transFlat.expanded());

        QScopedPointer<Cube, TemporaryCubeDeleter> translated(new Cube(transFlat.expanded(), "r"));

        flatcube.swap(translated);

        CubeAttributeInput att;

        p.SetInputCube(transFlat.expanded(),att);

      }  //Finished setting flatfield file for ONC-T

    }

    Cube *ocube  = p.SetOutputCube("TO");

    QString calfile = loadCalibrationVariables(ui.GetAsString("CONFIG"));

    g_timeRatio = g_Tvct/(g_texp + g_Tvct);

    QString g_units = "DN";

    // if ( "radiance" == g_iofCorrection.toLower() ) {

    //   // Units of RADIANCE

    calibrationLog.addKeyword(bnae);  

    calibrationLog.addKeyword(PvlKeyword("Bias_AETemp", toString(g_AEtemperature, 16)));

    switch (g_instrument) {

      case InstrumentType::ONCT:

         calibrationLog.addKeyword(PvlKeyword("Bias_CCDTemp", toString(g_CCD_T_temperature, 16)));

         calibrationLog.addKeyword(PvlKeyword("Bias_ECTTemp", toString(g_ECT_T_temperature, 16)));

    darkCurrent.addValue(toString(g_darkCurrent, 16));

    calibrationLog.addKeyword(darkCurrent);

    // Write Calibration group to output file

    ocube->putGroup(calibrationLog);

    Application::Log(calibrationLog);

  //configFile.clear();

    p.EndProcess();

  }

}