PROG: Cleaned up the formatting for amicacal.cpp/AmicaCalUtilities.h. References #3880,#4222.

#ifndef AmicaCalUtils_h

#define AmicaCalUtils_h

#include <cmath>

#include <string>

#include <vector>

#include "CSVReader.h"

#include "IException.h"

#include "IString.h"

#include "FileName.h"

#include "LineManager.h"

#include "NaifStatus.h"

#include "IString.h"

#include "Pvl.h"

#include "PvlGroup.h"

#include "Spice.h"

// OpenCV stuff

#include <opencv2/opencv.hpp>

/**

 * @author 2016-04-04 Tyler Wilson

 *

 *

 * @author 2016-07-10 Tyler Wilson

 * @internal

 *    @history

 */

using namespace cv;

using namespace std;

namespace Isis {

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

 * bodies.

 */

static void loadNaifTiming() {

  static bool naifLoaded = false;

  if (!naifLoaded) {

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

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

 */

static bool sunDistanceAU(const QString &scStartTime,

                          const QString &target,

                          double &sunDist) {

/**

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

    }

 }

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;

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

  }

 * at pixel with coordinates (x,y) relative to the central pixel (with coordinates (0,0) ).

 *

 */

static double f_focused(double alpha,int binning,double x,double y) {

        double X = x*binning;

 * central pixel (at coordinates (0,0) ).

 *

 */

static double f_unfocused(double * A,double * sigma, int N,int binning,double x,double y)  {

  double X = binning*x;

  double Y = binning*y;

  double r = sqrt(X*X+Y*Y);

  double sum = 0;

  for (int i = 0; i < N; i ++)   {     

    sum += (A[i]/(sigma[i]*sqrt(2.0*pi_c() ) ) )*exp(-(r*r)/(2*sigma[i]*sigma[i]) );

  }

  return sum;

}

 * @param binning

 * @return @b double * A pointer to a [size x size] matrix of light distribution values.

 */

double * setPSFFilter(int size, double *A,double *sigma, double alpha,int N,int binning) {

  double * psfVals = new double[size*size];

  int i = 0;

  for(double y = -(size / 2) ; y <= (size / 2) ; y++) {

    for(double x = -(size / 2) ; x <= (size / 2) ; x++) {

       if (x == 0 && y ==0) {

         psfVals[i] = 0;

         i++;

       }

       else {

         psfVals[i]=f_unfocused(A,sigma,N,binning,x,y) +f_focused(alpha,binning,x,y);

         //psfVals[i] = f_focused(alpha,binning,x,y);

         i++;

       }      

      //i++;

    }

  }

// $Id: amicacal.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 <vector>

#include <QFile>

#include <QString>

#include <QScopedPointer>

#include <QString>

#include <QTemporaryFile>

#include <QVector>

#include "iTime.h"

#include "LineManager.h"

#include "Pixel.h"

#include "ProcessByBoxcar.h"

#include "ProcessByBrick.h"

#include "ProcessByLine.h"

#include "ProcessBySample.h"

#include "ProcessByBrick.h"

#include "ProcessByBoxcar.h"

#include "ProgramLauncher.h"

#include "Pvl.h"

#include "PvlGroup.h"

//For subimage and binning mapping

static AlphaCube *alpha(0);

QString g_filter = "";

static QString g_target ="";

static int g_HayabusaNaifCode = -130;

static double g_temp(0);

static double g_darkCurrent(0);

//Smear calculation variables

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

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

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

  // Basic assurances...

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

    throw IException(IException::User,

  }

  try {

  g_temp = inst["CcdTemperature"] ;

  }

  catch(IException &e) {

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

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

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

  }

//tjw

  QString startTime = inst["SpacecraftClockStartCount"];

  g_startTime = startTime;

  binning = inst["Binning"];

  int lastline = inst["LastLine"];

  int lastsample = inst["LastSample"];

  nsubImages = archive["SubImageCount"];  // If > 1, some correction is 

                                          // not needed/performed

  nsubImages = archive["SubImageCount"];  /**<  If > 1, on-board smear-correction has been

                                           performed. */

  // I/F values

  QString target = inst["TargetName"];

  g_target = target;

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

    g_iof = pi_c()*(g_solarDist*g_solarDist)*(g_v_standard*g_iofScale);

  }

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

  //  occurred

  QScopedPointer<Cube, TemporaryCubeDeleter> flatcube;

  FileName flatfile= DetermineFlatFieldFile(g_filter,g_nullPolarizedPixels);

  QString reducedFlat(flatfile.expanded());

//Image is not cropped

if (startline ==0 && startsample == 0){

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

else {

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

      FileName transFlat =

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

  nl = icube->lineCount();

  nb = icube->bandCount();

  loadCalibrationVariables();

  g_timeRatio = g_Tvct/(g_texp+g_Tvct);

  g_darkCurrent = g_d0*exp(g_d1*g_temp);

  // Calibrate!

  ProcessByBoxcar pDiffusionModel;

  QScopedPointer<Cube, TemporaryCubeDeleter> diffusionModel;

  QString kernel_sz = QString::number(g_size);

  //QTemporaryFile psfModel("$TEMPORARY/psfModel.cub");

  QFile psfModel("psfModel"+QString::number(g_size)+".cub");

  //QFile("psfModel.cub");

  //psfModel.setAutoRemove(false);

  pDiffusionModel.SetInputCube(fname,attInput);

  pDiffusionModel.SetOutputCube(psfModel.fileName(),attOutput,ns,nl,nb);

  g_psfFilter = setPSFFilter(g_size, g_A,g_sigma, g_alpha,g_N,binning);

  try {

       pDiffusionModel.StartProcess(psfCorrectionBoxcar);  //Determine the diffusion model.

    }

    catch(IException &ie){

    throw IException(ie, IException::Programmer,

                     "Calculating the diffusion model failed!", _FILEINFO_);

    }

    pDiffusionModel.EndProcess();

    //Apply the PSF correction

      throw IException(ie, IException::Programmer,

                       "Applying the PSF correction failed!", _FILEINFO_);

    }    

   pPSFCorrection.EndProcess();

#endif

  // Log calibration activity

  PvlGroup calibrationLog("RadiometricCalibration");

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

/**

 * @brief Determine name of flat field file to apply

 * 

 * @author 2016-03-30 Kris Becker

 * 

 * @param filter  Name of AMICA filter

 * 

 * @return FileName Path and name of flat file file

 */

FileName DetermineFlatFieldFile(const QString &filter, const bool nullPolarPix) {

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

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

  if (nullPolarPix) {

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

  }

  else {

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

  }

  FileName final(fileName);

  //tjw:  So was this part

  //final = final.highestVersion();

  return final;

}

 * @param result  The light diffusion estimate at the central pixel of the boxcar of all the

 * surrounding pixels.

 */

void psfCorrectionBoxcar(Buffer &in, double &result) {

    result = 0;

    Statistics stats;

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

      if(!IsSpecial(in[i])) {

        stats.AddData(in[i]*g_psfFilter[i]);       

        //result += in[i]*g_psfFilter[i];

      }

    }

    result = stats.Sum();

}

 * @parm  in[1]  The radiometrically corrected cube (without PSF correction).

 * @param out    The radiometrically corrected cube after PSF correction has been applied.

 */

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

  Buffer& nopsf    = *in[1];

  Buffer& psfVals =  *in[0];

  Buffer& imageOut  = *out[0];

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

    if (!IsSpecial(psfVals[i])) {

    int samp(hotPixels[i][0].toInt());

    int line (hotPixels[i][1].toInt());

    hotPixelVector.append( Pixel(alpha->BetaSample(samp),alpha->BetaLine(line),1,0));

  }

  //Load linearity variables

  g_Gamma = Linearity["Gamma"];

  g_Gamma = 1.0-g_Gamma;

  //Load Smear Removal Variables

  g_Tvct = Smear["Tvct"];

  //Load DarkCurrent variables

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

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

 * @param in   Raw image and flat field

 * @param out  Radometrically corrected image

 */

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

  Buffer& imageIn   = *in[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;

  }

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

  double t1 = g_timeRatio/imageIn.size();

  }

  //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 (Isis::IsSpecial(imageOut[i])) {

      imageOut[i] = imageIn[i];

      }

    }

    // 2) LINEARITY Correction - always done

    if (imageOut[i] != Isis::Null) {

      imageOut[i] = pow(imageOut[i],g_Gamma) +g_L0*imageOut[i]*exp(g_L1*imageOut[i]);

    }

    // 3) DARK Current - Currently negligible and removed

#if 0

#endif

    // 4) HOT Pixel Removal

      bool hot = false;

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

        if ((hotPixelVector[j].sample() == currentSample) && (hotPixelVector[j].line() == i)) {

          imageOut[i] = Isis::Null;

          hot = true;

        }

      }

      if (hot == true)

        continue;

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

    //    accounted for in computation of c1 before loop.

    if ( imageOut[i] != Isis::Null && nsubImages <= 1 ) {

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

      }

    // 6) FLATFIELD correction

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

      if (Isis::IsSpecial(flatField[i])) {

        imageOut[i] = Isis::Null;

        continue;

          imageOut[i] /= flatField[i];

    }

    // 7) I/F Conversion

    if (imageOut[i] != Isis::Null && g_iofCorrection == true) {

          imageOut[i] *= g_iof;

       }