Added hyb2onccal calibration program.

#ifndef Hyb2OncCalUtils_h

#define Hyb2OncCalUtils_h

#include <cmath>

#include <string>

#include <vector>

#include <QString>

#include "CSVReader.h"

#include "IException.h"

#include "FileName.h"

#include "LineManager.h"

#include "NaifStatus.h"

#include "IString.h"

#include "Pvl.h"

#include "PvlGroup.h"

#include "Spice.h"

// OpenCV libraries

#include <opencv2/opencv.hpp>

/**

 * @author 2016-04-04 Tyler Wilson

 *

 *

 */

using namespace cv;

using namespace std;

namespace Isis {

/**

 * @brief Load required NAIF kernels required for timing needs.

 *

 * This method maintains the loading of kernels for HAYABUSA timing and

 * planetary body ephemerides to support time and relative positions of planet

 * bodies.

 */

/* Helper function for sunDistanceAu, don't need this until we have radiance calibration

   parameters for Hayabusa2 ONC-T filters to calculate radiance and I/F

static void loadNaifTiming() {

  static bool naifLoaded = false;

  if (!naifLoaded) {

//  Load the NAIF kernels to determine timing data

    Isis::FileName leapseconds("$base/kernels/lsk/naif????.tls");

    leapseconds = leapseconds.highestVersion();

    Isis::FileName sclk("$hayabusa2/kernels/sclk/hyb2_20141203-20161231_v01.tsc");    

    Isis::FileName pck1("$hayabusa2/kernels/tspk/de430.bsp");

    Isis::FileName pck2("$hayabusa2/kernels/tspk/jup329.bsp");

    Isis::FileName pck3("$hayabusa2/kernels/tspk/sat375.bsp");

    Isis::FileName pck4("$hayabusa2/kernels/spk/hyb2_20141203-20141214_0001m_final_ver1.oem.bsp");

    Isis::FileName pck5("$hayabusa2/kernels/spk/hyb2_20141203-20151231_0001h_final_ver1.oem.bsp");

    Isis::FileName pck6("$hayabusa2/kernels/spk/hyb2_20151123-20151213_0001m_final_ver1.oem.bsp");

//  Load the kernels

    QString leapsecondsName(leapseconds.expanded());

    QString sclkName(sclk.expanded());

    QString pckName1(pck1.expanded());

    QString pckName2(pck2.expanded());

    QString pckName3(pck3.expanded());

    QString pckName4(pck4.expanded());

    QString pckName5(pck5.expanded());

    QString pckName6(pck6.expanded());

    furnsh_c(leapsecondsName.toLatin1().data());

    furnsh_c(sclkName.toLatin1().data());

    furnsh_c(pckName1.toLatin1().data());

    furnsh_c(pckName2.toLatin1().data());

    furnsh_c(pckName3.toLatin1().data());

    furnsh_c(pckName4.toLatin1().data());

    furnsh_c(pckName5.toLatin1().data());

    furnsh_c(pckName6.toLatin1().data());

//  Ensure it is loaded only once

    naifLoaded = true;

  }

  return;

}

*/

/**

 * @brief Computes the distance from the Sun to the observed body.

 *

 * This method requires the appropriate NAIK kernels to be loaded that

 * provides instrument time support, leap seconds and planet body ephemeris.

 *  

 * @return @b double Distance in AU between Sun and observed body.

 */

/* commented out until we have radiance values (RAD/IOF group in calibration trn) for Hayabusa2.

 static bool sunDistanceAU(const QString &scStartTime,

                          const QString &target,

                          double &sunDist) {

  //  Ensure NAIF kernels are loaded

  loadNaifTiming();

  sunDist = 1.0;

  //  Determine if the target is a valid NAIF target

  SpiceInt tcode;

  SpiceBoolean found;

  bodn2c_c(target.toLatin1().data(), &tcode, &found);

  if (!found) return (false);

  //  Convert starttime to et

  double obsStartTime;

  scs2e_c(-37, scStartTime.toLatin1().data(), &obsStartTime);

  //  Get the vector from target to sun and determine its length

  double sunv[3];

  double lt;

  spkpos_c(target.toLatin1().data(), obsStartTime, "J2000", "LT+S", "sun",

                  sunv, &lt);

  NaifStatus::CheckErrors();

  double sunkm = vnorm_c(sunv);

  //  Return in AU units

  sunDist = sunkm / 1.49597870691E8;

  //cout << "sunDist = " << sunDist << endl;

  return (true);

}

*/

/**

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

 *

 * @author 2016-04-19 Tyler Wilson

 *

 * @param icube A pointer to the input cube

 *

 * @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++) {

    linereader.SetLine(line+1);

    icube->read(linereader);

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

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

    }

 }

return matrix;

}

/**

 * @brief Translates an OpenMat object to an ISIS::Cube with one band

 *

 * @author 2016-04-19 Tyler Wilson

 *

 * @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) {

  int nlines = matrix->rows;

  int nsamples = matrix->cols;

  CubeAttributeOutput set;

  set.setPixelType(Real);

  Cube ocube;

  ocube.setDimensions(nsamples,nlines,1);

  ocube.create(cubeName,set);

  LineManager linewriter(ocube);

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

    linewriter.SetLine(line+1);

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

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

    }

    ocube.write(linewriter);

  }

}

/**

 * @brief Translates/scales a cube using Bilinear Interpolation

 *

 * @author 2016-04-19 Tyler Wilson

 *

 * @param matrix A pointer to the OpenMat object

 *

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

 *

 */

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

  Mat * originalMat = isis2mat(flatField);

  int scale = transform[0];

  int startsample = transform[1];

  int startline = transform[2];

  int lastsample = transform[3];

  int lastline = transform[4];

  int width  = (lastsample-startsample);

  int height = (lastline-startline);

  Size sz(flatField->lineCount()/scale,flatField->sampleCount()/scale);

  Mat * resizedMatrix = new Mat();

  Mat temp = *originalMat;

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

  if (scale ==1) {   

    mat2isis(&originalCropped,fname);

  }

  else {

    //Bilinear interpolation

    resize(originalCropped,*resizedMatrix,sz,INTER_LINEAR);

    mat2isis(resizedMatrix,fname);

  }

}

}

#endif

ifeq ($(ISISROOT), $(BLANK))

.SILENT:

error:

	echo "Please set ISISROOT";

else

	include $(ISISROOT)/make/isismake.apps

endif

# opencv includes path and libs that are not in

# isismake.os, but we don't want all apps to 

# link against the large opencv libraries

# so they are added here:

#

ALLLIBS += $(OPENCVLIBS)

// $Id: hyb2onccal.cpp 6045 2015-02-07 02:06:59Z moses@GS.DOI.NET $

#include "Isis.h"

#include <vector>

#include <algorithm>

#include <memory>

#include <cstdio>

#include <cmath>

#include <QFile>

#include <QString>

#include <QScopedPointer>

#include <QTemporaryFile>

#include <QVector>

#include "AlphaCube.h"

#include "Buffer.h"

#include "FileName.h"

#include "Hyb2OncCalUtils.h"

#include "IException.h"

#include "iTime.h"

#include "LineManager.h"

#include "Pixel.h"

#include "ProcessByLine.h"

#include "ProcessBySample.h"

#include "ProcessByBrick.h"

#include "ProcessByBoxcar.h"

#include "ProgramLauncher.h"

#include "Pvl.h"

#include "PvlGroup.h"

#include "Spice.h"

#include "Statistics.h"

#include "TextFile.h"

using namespace Isis;

using namespace std;

// Calibration support routines

FileName DetermineFlatFieldFile(const QString &filter);

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

QString loadCalibrationVariables(const QString &config);

// 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() );

    }

  }

};

//For subimage and binning mapping

static AlphaCube *alpha(0);

QString g_filter = "";

static QString g_target ="";

static const int g_Hayabusa2NaifCode = -130;

static Pvl g_configFile;

//Bias calculation variables

static double g_b0(0);

static double g_b1(0);

static double g_b2(0);

static double g_bias(0);

static double g_AEtemperature(0);

static QString g_startTime;

//Dark Current variables

static double g_d0(0);

static double g_d1(0);

static double g_temp(0);

static double g_darkCurrent(0);

// TODO: we do not have the readout time (transfer period) for Hayabusa2 ONC.

//Smear calculation variables

// static double g_Tvct(0);       //!< Vertical charge-transfer period (in seconds).

static double g_texp(1);       //!< Exposure time.

// static double g_timeRatio(1.0);

// Calibration parameters

static int nsubImages(0);      //!< Number of sub images

static int binning(1);         //!< The number of samples/lines which are binned

static double g_compfactor(1.0);  // Default if OutputMode = LOSS-LESS; 16.0 for LOSSY

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

//  I/F variables

static double g_solarDist(1.0);  /**< Distance from the Sun to the target body

(used to calculate g_iof) */

static double g_iof(1.0);        //!< I/F conversion value

static double g_iofScale(1.0);

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

// TODO: we do not have this conversion factor for Hayabusa 2 ONC.

static double g_v_standard(1.0);

// static double g_v_standard(3.42E-3);//!< Base conversion for all filters (Tbl. 9)

void IsisMain() {

  UserInterface& ui = Application::GetUserInterface();

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

  const QString hyb2cal_program = "hyb2onccal";

  const QString hyb2cal_version = "1.0";

  const QString hyb2cal_revision = "$Revision$";

  QString hyb2cal_runtime = Application::DateTime();

  ProcessBySample p;

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

  // Basic assurances...

  if (icube->bandCount() != 1) {

    throw IException(IException::User,

      "ONC images may only contain one band", _FILEINFO_);

    }

  PvlGroup &inst = icube->group("Instrument");

  PvlGroup &bandbin = icube->group("BandBin");

  QString filter = bandbin["FilterName"];

  g_filter = filter;

  //Set up binning and image subarea mapping

  binning = inst["Binning"];

  int startLine = inst["SelectedImageAreaY1"];

  int startSample = inst["SelectedImageAreaX1"];

  int lastLine = inst["SelectedImageAreaY4"];

  int lastSample = inst["SelectedImageAreaX4"];

  AlphaCube myAlpha(1024,1024,icube->sampleCount(), icube->lineCount(),

  startSample,startLine,lastSample,lastLine);

  alpha = &myAlpha;

  try {

    g_texp = inst["ExposureDuration"] ;

  }

  catch(IException &e) {

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

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

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

  }

  QString instrumentId = inst["InstrumentId"];

  QString oncCCDTemperature = "ONC";

  try {

    oncCCDTemperature += instrumentId.section('-',1,1) + "CCDTemperature";

    g_temp = inst[oncCCDTemperature];

  }

  catch(IException &e) {

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

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

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

  }

  QString startTime = inst["SpacecraftClockStartCount"];

  g_startTime = startTime;

  nsubImages = inst["SubImageCount"];  // If > 1, some correction is not needed/performed

  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";

  if (instrumentId == "ONC-T") {

    QScopedPointer<Cube, TemporaryCubeDeleter> flatcube;

    flatfile = DetermineFlatFieldFile(g_filter);

    QString reducedFlat(flatfile.expanded());

    // Image is not cropped

    if (startLine == 1 && startSample == 1) {

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

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

      if (binning > 1) {

        QString scale(toString(binning));

        FileName newflat = FileName::createTempFile("$TEMPORARY/" +

        flatfile.baseName() + "_reduced.cub");

        reducedFlat = newflat.expanded();

        QString parameters = "FROM=" + flatfile.expanded() +

        " TO="   + newflat.expanded() +

        " MODE=SCALE" +

        " LSCALE=" + scale +

        " SSCALE=" + scale;

        try {

          ProgramLauncher::RunIsisProgram("reduce", parameters);

        }

        catch (IException& ie) {

          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);

    }

    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() );

      int transform[5] = {binning,startSample,startLine,lastSample,lastLine};

      // 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;

      std::cout<<" before second setinputcube()...\n"<<std::endl;

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

      std::cout<<" after second setinputcube()...\n"<<std::endl;

    }

  }

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

  QString fname = ocube->fileName();

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

  // TODO: we do not have g_Tvct (readout time for the ccd). Once we know this value for Hayabusa2

  // ONC and we know the smear model, we can correctly account for smear.

  //g_timeRatio = g_Tvct/(g_texp + g_Tvct);

  g_darkCurrent = g_d0 * exp(g_d1 *g_temp);

  std::cout << "\ndark current: " << g_darkCurrent << std::endl;

  g_iof = 1.0;  // Units of DN

  QString g_units = "DN";

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

  //   // Units of RADIANCE

  //   g_iof = g_iof * g_iofScale;

  //   g_units = "W / (m**2 micrometer sr)";

  // }

  //

  // if ( !sunDistanceAU(startTime, target, g_solarDist) ) {

  //    throw IException(IException::Programmer, "Cannot calculate distance to sun!",

  //                      _FILEINFO_);

  // }

  //

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

  //   // Units of I/F

  //   // TODO: Note, we do not have a correct g_iofScale (== 1 right now). This was provided for

  //   // Hayabusa AMICA's v-band and all other bands were normalized according to this value. We do

  //   // not have this value for Hayabusa2 ONC. Need to correct this value.

  //   g_iof = pi_c() * (g_solarDist * g_solarDist) *

  //           g_iofScale / g_solarFlux / g_texp;

  //   g_units = "I over F";

  // }

  // Calibrate!

  try {

    p.Progress()->SetText("Calibrating Hayabusa2 Cube");

    p.StartProcess(Calibrate);

  }

  catch (IException &ie) {

    throw IException(ie, IException::Programmer,

      "Radiometric calibration failed!", _FILEINFO_);

  }

  // Log calibration activity performed so far

  PvlGroup calibrationLog("RadiometricCalibration");

  calibrationLog.addKeyword(PvlKeyword("SoftwareName", hyb2cal_program));

  calibrationLog.addKeyword(PvlKeyword("SoftwareVersion", hyb2cal_version));

  calibrationLog.addKeyword(PvlKeyword("ProcessDate", hyb2cal_runtime));

  calibrationLog.addKeyword(PvlKeyword("CalibrationFile", calfile));

  calibrationLog.addKeyword(PvlKeyword("FlatFieldFile", flatfile.originalPath()

  + "/" + flatfile.name()));

  calibrationLog.addKeyword(PvlKeyword("CompressionFactor", toString(g_compfactor, 2)));

  // Parameters

  PvlKeyword key("Bias_Bn");

  key.addValue(toString(g_b0, 8));

  key.addValue(toString(g_b1, 8));

  key.addValue(toString(g_b2, 8));

  calibrationLog.addKeyword(key);

  calibrationLog.addKeyword(PvlKeyword("Bias", toString(g_bias, 16), "DN"));

  // calibrationLog.addKeyword(PvlKeyword("Smear_Tvct", toString(g_Tvct, 16)));

  calibrationLog.addKeyword(PvlKeyword("CalibrationUnits", g_iofCorrection));

  calibrationLog.addKeyword(PvlKeyword("RadianceStandard", toString(g_v_standard, 16)));

  calibrationLog.addKeyword(PvlKeyword("RadianceScaleFactor", toString(g_iofScale, 16)));

  calibrationLog.addKeyword(PvlKeyword("SolarDistance", toString(g_solarDist, 16), "AU"));

  calibrationLog.addKeyword(PvlKeyword("SolarFlux", toString(g_solarFlux, 16)));

  calibrationLog.addKeyword(PvlKeyword("IOFFactor", toString(g_iof, 16)));

  calibrationLog.addKeyword(PvlKeyword("Units", g_units));

  // Write Calibration group to output file

  ocube->putGroup(calibrationLog);

  Application::Log(calibrationLog);

  //configFile.clear();

  p.EndProcess();

}

/**

* @brief Determine name of flat field file to apply

* @author 2016-03-30 Kris Becker

* @param filter  Name of ONC filter

* @return FileName Path and name of flat file file

* @internal

*   @history 2017-07-27 Ian Humphrey & Kaj Williams - Adapted from amicacal.

*/

FileName DetermineFlatFieldFile(const QString &filter) {

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

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

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

  FileName final(fileName);

  return final;

}

/**

* @brief Loads the calibration variables into the program.

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

*/

QString loadCalibrationVariables(const QString &config)  {

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

  FileName calibFile(config);

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

  // Pvl configFile;

  g_configFile.read(calibFile.expanded());

  // Load the groups

  PvlGroup &Bias = g_configFile.findGroup("Bias");

  PvlGroup &DarkCurrent = g_configFile.findGroup("DarkCurrent");

  // PvlGroup &Smear = g_configFile.findGroup("SmearRemoval");

  PvlGroup &solar = g_configFile.findGroup("SOLARFLUX");

  // PvlGroup &iof = g_configFile.findGroup("RAD");

  // Load Smear Removal Variables

  // g_Tvct = Smear["Tvct"];

  // Load DarkCurrent variables

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

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

  // Load Bias variables

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

  g_b1 = Bias["B"][1].toDouble();

  g_b2 = Bias["B"][2].toDouble();

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

  g_bias = g_b0 + g_b1 * g_AEtemperature + g_b2 * (g_AEtemperature * g_AEtemperature);

  std::cout << "bias: " << g_bias << std::endl;

  // Load the Solar Flux for the specific filter

  g_solarFlux = solar[g_filter.toLower()];

  // radiance = g_v_standard * g_iofScale

  // iof      = radiance * pi *dist_au^2

  // g_iofScale   = iof[g_filter];

  return ( calibFile.original() );

}

/**

* @brief Apply radiometric correction to each line of an AMICA image.

* @author 2016-03-30 Kris Becker

* @param in   Raw image and flat field

* @param out  Radometrically corrected image

* @internal

*   @history 2017-07-2017 Ian Humphrey & Kaj Williams - Adapted from amicacal.

*/

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

  Buffer& imageIn   = *in[0];

  Buffer& flatField = *in[1];

  Buffer& imageOut  = *out[0];

  int pixelsToNull = 0;

  // Note that is isn't currently tested, as we do not have a test with a hayabusa2 image that

  // has been on-board cropped.

  int currentSample = imageIn.Sample();

  int alphaSample = alpha->AlphaSample(currentSample);

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

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

      imageOut[i] = Isis::Null;

    }

    return;

  }

  // TODO: Once smear model and readout time are known, we can add smear correction.

  // Compute smear component here as its a constant for the entire sample

  // double t1 = g_timeRatio/imageIn.size();

  // double b = binning;

  // double c1(1.0);  //default if no binning

  //

  // if (binning > 1) {

  //   c1 = 1.0/(1.0 + t1*((b -1.0)/(2.0*b) ) );

  // }

  //

  // double smear = 0;

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

  //   if ( !IsSpecial(imageIn[j]) ) {

  //     smear += t1 * ( (imageIn[j] * g_compfactor) - g_bias);

  //   }

  // }

  // Iterate over the line space

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

    imageOut[i] = imageIn[i];

    // Check for special pixel in input image and pass through

    if ( IsSpecial(imageOut[i]) ) {

      imageOut[i] = imageIn[i];

      continue;

    }

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

    imageOut[i] *= g_compfactor;

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

    // TODO: find image with > 1 subimages

    if ( nsubImages <= 1 ) {

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

        imageOut[i] = Null;

        continue;

      }

      else {

        imageOut[i] = imageOut[i] - g_bias;

      }

    }

    // 2) DARK Current

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

    // READOUT Smear Removal - Not needed if on-board corrected.  Binning is

    //    accounted for in computation of c1 before loop.

    // if (nsubImages <= 1) {

    //  imageOut[i] = c1*(imageOut[i] - smear);

    // }

    // 3) FLATFIELD correction

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

    // If we have only one input cube, that means that we do not have a flat-field (W1/W2).

    if (in.size() == 2) {

      // Note that our current flat-fields to not have special pixel values.

      if ( IsSpecial(flatField[i]) ) {

        imageOut[i] = Isis::Null;

        continue;

      }

      else {

        if (flatField[i] != 0) {

          imageOut[i] /= flatField[i];

        }

      }

    }

    // TODO: once the radiance values are known for each band, we can correctly compute I/F.

    // For now, g_iof is 1, so output will be in DNs.

    // 7) I/F or Radiance Conversion (or g_iof might = 1, in which case the output will be in DNs)

    imageOut[i] *= g_iof;