Added lromakeflat application from LRO mission team.

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

.SILENT:

error:

	echo "Please set ISISROOT";

else

	include $(ISISROOT)/make/isismake.apps

endif

 No newline at end of file

/**

 * @Brief This application creates three flatfield (Sensitivity Non-Uniformity Matrix) cubes used for calibration

 *

 * This application creates three flatfield (Sensitivity Non-Uniformity Matrix) cubes used for calibration.

 * The cubes consist of median, mean, and standard deviation values per pixel. Process varies for the three 

 * cameras this application can be used for but general pixel stacking column approach is the same. The 

 * three camera types are line-scan, push-frame, and framing. Invalid pixel values are changed to Isis::Null. 

 *

 * The application uses a two step process.

 *

 * Step #1 - STACKING PIXELS INTO PIXEL COLUMNS

 *

 * The first part of the process is to stack cube pixels at their respective pixel locations.

 *

 * These stacked pixels are processed to produce the median, mean, and standard deviation values for the 

 * pixel location, and to toss out invalid pixel values.  Normalization can also be done at this time if 

 * images were not normalized prior. The mean is used to normalize images. 

 *

 * Step #2 - PROCESSING STACKED PIXEL COLUMNS AND TRANSFERRING TO SNU

 *

 * Once pixels have been stacked into pixel columns mean and median averages, and standard devation for that 

 * pixel location is saved to the SNU matrix (flatfield cubes).

 *

 * Mean

 *  -The summed pixels are divided by the number of pixels.

 * Standard Deviation 

 *  -The mean average is used to calculate the standard deviation. 

 * Median

 *  -The pixel column is sorted and the median is obtained.

 *

 * @author Victor Silva - 2016-08-17

 *

 * @internal

 *   @history 2016-08-17 Victor Silva - New app for lroc adapted from makeflat app

 *

 */

#include "Isis.h"

#include <bits/algorithmfwd.h>

#include <cmath>

#include <iostream>

#include <limits>

#include <new>

#include <string>

#include <stdlib.h>

#include <vector>

#include <QString>

#include "FileList.h"

#include "Histogram.h"

#include "IException.h"

#include "LineManager.h"

#include "PvlGroup.h"

using namespace std;

using namespace Isis;

// Container for two-dimensional matrix of doubles

typedef vector < vector<double> > Matrix2d;

// This method returns stats for a vector of doubles and excludes special pixels

void getVectorStats( vector<double> (&inVec), double (&mean), double (&median), double (&stdev) );

// Get normalization values for frames in list of cubes if so desired

void getCubeListNormalization( Matrix2d (&matrix), int (cubeWidth), int (frameHeight), long (frameLineCount), long (iLineCount) );

// Logs excluded files and frames

void exclude( int listIndex, string reason );

// Logs excluded files and frames

void exclude( int listIndex, int frame, string reason );

// Global variables

int static g_sampleCount;

int static g_numStdevs;

QString static g_msg;

Progress static g_prog;

FileList static g_list;

// This will keep track of excluded cubes from the cube list. The int is the index in the list (zero-based)

map < int, string > g_excludedCubes;

// Holds excluded frames. <cube index in list, frame index in cube> (zero-based)

map < pair<int, int>, string> g_excludedFrames;

// Holds exclusion details used for logging

vector < PvlObject > g_excludedDetails;

/**

   * @brief Main method

   *

   * Computes flatfield image for linescan, pushbroom, and framing

   * Cameras

   * 

   * @throws IException::User "Only single band images accepted"

   * @throws IException::User "User selected lines value exceeds number of lines in cube"

   * @throws IException::Programmer "Could not write to output cube"

   *

   * @author 2016-08-17 Victor Silva

   *

   * @internal

   *   @history 2016-08-17 Victor Silva - New app for lroc adapted from makeflat app

   *                         

   */

void IsisMain() {

  // Set attributes

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

  g_list = FileName(ui.GetFileName("FROMLIST"));

  g_numStdevs = abs(ui.GetDouble("STDEVTOL"));

  // Input cube attributes

  bool normalize = false;

  long iLineCount, oLineCount; // sample count is the same for both input/output cubes

  int cubeWidthPixels = 0, frameHeightLines = 0, framesPerCube, frameLineCount;

  // Setup Sensitivity Non-uniformity Matrix cubes (flatfielding matrix) and line managers

  Cube *oStdevCube = NULL, *oMedianCube = NULL, *oMeanCube = NULL; 

  LineManager *oStdevlineMgr = NULL, *oMeanLineMgr = NULL, *oMedianLineMgr = NULL;

  // Get some info from first cube in list

  Cube firstCube;

  firstCube.open(g_list[0].toString());

  g_sampleCount = firstCube.sampleCount();

  iLineCount = firstCube.lineCount();

  if (firstCube.bandCount() != 1) {

    string err = "Only single band images accepted. \n";

    throw IException(IException::User, err, _FILEINFO_);

  }

  firstCube.close();

  // to normalize or not

  if (ui.GetString("NORMALIZATION") == "YES") {

      normalize = true;

      frameHeightLines = abs(ui.GetInteger("HEIGHT_LINES"));

      cubeWidthPixels = abs(ui.GetInteger("WIDTH_PIXELS"));

  }

  // User selection of camera type will set addt'l attributes

  if (ui.GetString("IMAGETYPE") == "FRAMING") {

    frameLineCount = oLineCount = iLineCount;

  }

  else if (ui.GetString("IMAGETYPE") == "LINESCAN") {

    oLineCount = 1;

    frameLineCount = abs(ui.GetInteger("NUMLINES"));

    if (iLineCount != frameLineCount) {

      string err = "User selected lines value (" + IString(frameLineCount) + 

         ") exceeds number of lines in cube (" + IString(iLineCount) + ". \n";

      throw IException(IException::User, err, _FILEINFO_);

    }

  }

  else {

    oLineCount = frameLineCount = abs(ui.GetInteger("FRAMELETHEIGHT"));

  }

  // Number of frames in each cube

  framesPerCube = (iLineCount / frameLineCount); //integer division

  // 3d matrix of type pair to hold mean and stdev for normalization

  Matrix2d normMatrix;

  normMatrix.resize(g_list.size(), vector<double> (framesPerCube, 1.0));

  // Get normalization values from list of cubes if user selects normalization

  if (normalize) {

    getCubeListNormalization(normMatrix, cubeWidthPixels, frameHeightLines, frameLineCount, iLineCount);

  }

  // Create empty snu-flatfield output cubes and line managers

  oStdevCube = new Cube();

  oMedianCube = new Cube();

  oMeanCube = new Cube();

  // Set cube dimensions

  oStdevCube->setDimensions(g_sampleCount, oLineCount, 1);

  oMedianCube->setDimensions(g_sampleCount, oLineCount, 1);

  oMeanCube->setDimensions(g_sampleCount, oLineCount, 1);

  // Set name of cubes

  oStdevCube->create(FileName(ui.GetFileName("TO")).expanded() + ".stdev.cub");

  oMedianCube->create(FileName(ui.GetFileName("TO")).expanded() + ".median.cub");

  oMeanCube->create(FileName(ui.GetFileName("TO")).expanded() + ".mean.cub");

  // create new line manager

  oStdevlineMgr = new LineManager(*oStdevCube);

  oMedianLineMgr = new LineManager(*oMedianCube);

  oMeanLineMgr = new LineManager(*oMeanCube);

  // Set status msg

  g_msg = "Stacking pixels into columns and processing them.";

  g_prog.SetText(g_msg);

  g_prog.SetMaximumSteps(oLineCount + 1);

  // Stack pixels at each pixel location into pixel columns

  for (long oLine = 1; oLine <= oLineCount; oLine++) {

    g_prog.CheckStatus();

    // 2d matrix to hold pixel columns

    Matrix2d pixelMatrix;

    pixelMatrix.resize(g_sampleCount, vector<double>(g_list.size() * framesPerCube, Isis::Null));

    for (int listIndex = 0; listIndex < g_list.size(); listIndex++) {

      Cube tmp2;

      tmp2.open(g_list[listIndex].toString());

      // Only run for cubes with one band

      if (tmp2.bandCount() != 1) {

        string err = "Warning: This cube has too many bands(" + IString(tmp2.bandCount()) + 

          " and will be excluded). Only single band images accepted. \n";

        cerr << err << endl;

        exclude(listIndex, err);

        // Go directly to next cube. Do not stop at GO. Do not collect $200.

        continue;

      }

      // Reset frame for every cube in list

      int frame = 0;

      // Cube line has to match oLine as we are slicing the cubes one line at a time to stack pixels 

      for (long cubeLine = oLine; cubeLine <= iLineCount; cubeLine = cubeLine + frameLineCount) {

        LineManager cubeMgr(tmp2);

        cubeMgr.SetLine(cubeLine);

        tmp2.read(cubeMgr);

        double pixelVal;

        // Iterate through columns

        for (int column = 0; column < g_sampleCount; column++) {

          pixelVal = cubeMgr[column];

          int frameIndex = listIndex * framesPerCube + frame;

          // If frame avg is zero, will cause div/zero error. Set pixel to 0? nil?  

          (normMatrix[listIndex][frame])? pixelVal = (pixelVal/(normMatrix[listIndex][frame])):pixelVal = Isis::Null;

          pixelMatrix[column][frameIndex] = pixelVal;

        }

        // Next frame

        frame++;

      }

      tmp2.close();

    }

    // Process stacked pixel columns and write to snu-flatfield output cubes

    // Set ouput cubeLine manager line to cube line

    oStdevlineMgr->SetLine(oLine);

    oMedianLineMgr->SetLine(oLine);

    oMeanLineMgr->SetLine(oLine);

    // Set values as we iterate through columns

    for (int column = 0; column < g_sampleCount; column++) {

      double columnMean, columnMedian, columnStdev;

      getVectorStats(pixelMatrix[column], columnMean, columnMedian, columnStdev);

      (*oStdevlineMgr)[column] = columnStdev;

      (*oMedianLineMgr)[column] = columnMedian;

      (*oMeanLineMgr)[column] = columnMean;

    }

    try {

      // Write the stats saved at the pixel location to the snu matrix cubes

      oStdevCube->write(*oStdevlineMgr);

      oMeanCube->write(*oMeanLineMgr);

      oMedianCube->write(*oMedianLineMgr);

    }

    catch (IException e) {

      //cout << "It's no use at this point. Just go. Leave me here. I'll only slow you down.";

      string err = "Could not write to output cube " + IString(FileName(ui.GetFileName("TO")).expanded()) + ".\n";

      throw IException(IException::Programmer, err, _FILEINFO_);

    }

    pixelMatrix.clear();

  }

  // Clean-up

  normMatrix.clear();

  if (oStdevlineMgr ){

    delete oStdevlineMgr;

    oStdevlineMgr = NULL;

  }

  if(oMeanLineMgr){

    delete oMeanLineMgr;

    oMeanLineMgr = NULL;

  }

  if(oMedianLineMgr){

    delete oMedianLineMgr;

    oMedianLineMgr = NULL;

  }

  if(oStdevCube) {

    oStdevCube->close();

    delete oStdevCube;

    oStdevCube = NULL;

  }

  if(oMeanCube) {

    oMeanCube->close();

    delete oMeanCube;

    oMeanCube = NULL;

  }

  if(oMedianCube) {

    oMedianCube->close();

    delete oMedianCube;

    oMedianCube = NULL;

  }

  // If user wanted an exclusion file created, create it pues.

  if(ui.WasEntered("EXCLUDE")) {

    Pvl excludeFile;

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

      excludeFile.addObject(g_excludedDetails[i]);

    }

    excludeFile.write(FileName(ui.GetFileName("EXCLUDE")).expanded());

  }

} // end of main

/**

 * @brief This function gets normalization values for a list of cubes

 *

 * It uses the user specified area of the image to be used for normalizing all

 * the values at each pixel location. The user enters a height and width

 * percentage of the frame.

 *

 * @param matrix  3D matrix that holds cube frame information such as stdev/normAvg

 * @param int cubeWidth size of cube in samples

 * @param int frameHeight size of frame in rows

 * @param long frameLineCount number of lines to make a frame

 * @param long iLineCount number of lines total in cube

 * 

 * @throws IException::User "This selection will yield less than 1 pixel (width). This is not enough to normalize."

 * @throws IException::User "This percentage will yield less than 1 line. This is not enough to normalize."

 *

 * @author 2016-08-17 Victor Silva

 *

 * @internal

 *   @history 2016-08-17 Victor Silva - New app for lroc adapted from makeflat app

 *   

 */

void getCubeListNormalization(Matrix2d &matrix, int cubeWidth, int frameHeight, long frameLineCount, long iLineCount) {

  for(int listIndex = 0; listIndex < g_list.size(); listIndex++) {

    // open cubes in list unless they have already been excluded

    if( g_excludedCubes.find(listIndex) != g_excludedCubes.end() ) {

      continue;

    }

    Cube tmp;

    tmp.open(g_list[listIndex].toString());

    PvlGroup normResults("NormalizationResults");

    g_msg = "Getting frame mean avg " + toString((int) listIndex + 1) + "/";

    g_msg += toString((int) g_list.size()) + " (" + g_list[listIndex].name() + ")";

    g_prog.SetText(g_msg);

    g_prog.SetMaximumSteps(g_list.size());

    g_prog.CheckStatus();

    // width:

    int startSample = 0, endSample = g_sampleCount;

    if (cubeWidth){

      if(cubeWidth <= g_sampleCount) {

        startSample = (g_sampleCount - cubeWidth) / 2;

        endSample = (g_sampleCount - 1) - startSample;

      }

      else {

        string err = "This selection will yield less than 1 pixel (width). This is not enough to normalize. \n";

        throw IException(IException::User, err, _FILEINFO_);

      }

    }

    // height:

    long startLine, endLine;

    int lineShift = 0;

    (frameHeight == 0) ? frameHeight = frameLineCount : lineShift = (frameLineCount - frameHeight)/ 2;

    for(long cubeLine = 1; cubeLine <= iLineCount; cubeLine = cubeLine + frameLineCount) {

      startLine = cubeLine + lineShift;

      endLine = startLine + frameHeight;

      if( (startLine < 0)||(endLine - startLine) < 0) {

        string err = "This percentage will yield less than 1 line. This is not enough to normalize. \n";

        throw IException(IException::User, err, _FILEINFO_);

      }

      Histogram hist = Histogram(tmp, 1, &g_prog, startSample, startLine, endSample, endLine, 0, true);

      int frame = cubeLine/frameLineCount;

      double normalizationAverage = hist.Average();

      matrix[listIndex][frame] = normalizationAverage;

      normResults += PvlKeyword("FileName", g_list[listIndex].toString());

      normResults += PvlKeyword("Frame", toString(frame));

      normResults += PvlKeyword("Frame_MeanAverage", toString(normalizationAverage));

      Application::Log(normResults);

    }

    tmp.close();

  }

}

/**

 * @brief This function returns the median, mean, and standard deviation for a vector 

 *        of valid pixels

 *

 * It will return Isis::Null if it doesn't find any valid pixels. The function builds 

 * a vector of valid pixels. It also calculates the mean and standard deviation of the 

 * vector. The function also sorts the vector and finds the median. If there are an even

 * number of valid pixels, it will average (mean) the two middle pixels and return

 * that average as the median. It does not return the median pixel but rather its

 * value or average value of the two middle pixels.

 *

 * @param inVec   Vector of pixel values (double)

 * @param mean    Value of vector mean (double)

 * @param median  Value of vector median (double)

 * @param stdev   Value of vector stdev (double)

 *

 * @author 2016-08-17 Victor Silva

 *

 * @internal

 *   @history 2016-08-17 Victor Silva - New app for lroc adapted from makeflat app

 */

void getVectorStats(vector<double> &inVec, double &mean, double &median, double &stdev) {

  // Set defaults

  int medianIndex = 0, v_size = 0;

  double sum = 0.0, devSum = 0.0;

  median = Isis::Null;

  mean = Isis::Null;

  stdev = Isis::Null;

  // Check for empty initial vector

  if(!(inVec.empty())) {

    vector <double> v1;

    v_size = int(inVec.size());

    // strip invalid pixels

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

      if (!(IsNullPixel(inVec[i]) || IsSpecial(inVec[i]))){

        v1.push_back(inVec[i]);

        sum += inVec[i];

      }

    }

    // check to see if we got any valid pixels

    if (!(v1.empty())) {

      v_size = int(v1.size());

      // First pass - toss out outliers

      double tempMean = sum / v_size;

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

        devSum += pow((v1[i] - tempMean), 2);

      }

      double tempStdev = sqrt(devSum / (v_size));

      // Reset for second pass

      devSum = 0;

      sum = 0;

      vector <double> v2;

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

        if( (abs(tempMean - v1[i])) <= (tempStdev * g_numStdevs) ){  

          v2.push_back(v1[i]);

          sum += v1[i];

        }

      }

      // Second pass - get stats

      if (!(v2.empty())) {

        v_size = int(v2.size());

        // mean

        mean = sum / v_size;

        // stdev process

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

          devSum += pow((v2[i] - mean), 2);

        }

        stdev= sqrt(devSum / (v_size));

        // median

        bool even = ((v_size % 2) == 0);

        sort(v2.begin(), v2.end());

        medianIndex = (v_size - 1) / 2;

        if (even) {

          median = double(((v2[medianIndex]) + (v2[medianIndex + 1])) / 2.0000);

        }

        else {

          median = double(v2[medianIndex]);

        }

        v2.clear();

      } 

      v1.clear();

    } 

  } 

}

/**

 * @brief This function excludes cubes and logs them

 * Adds excluded cube files to containers that are used when creating the

 * exclusion file. It also logs the exclusions to the application log window.

 *

 * @param listIndex  Location in list for a cube (zero-indexed)

 * @param reason     Exclusion reason

 *

 * @author 2016-08-17 Victor Silva

 *

 * @internal

 *   @history 2016-08-17 Victor Silva - New app for lroc adapted from makeflat app

 */

void exclude(int listIndex, string reason) {

  if(g_excludedCubes.find(listIndex) == g_excludedCubes.end()) {

    g_excludedCubes.insert(pair<int, string>(listIndex, reason));

    PvlObject exclusion("Excluded_Items");

    PvlGroup excludedFiles("Excluded_Files");

    excludedFiles += PvlKeyword("FileName", g_list[listIndex].toString());

    excludedFiles += PvlKeyword("Reason", QString::fromStdString(reason));

    Application::Log(excludedFiles);

    exclusion.addGroup(excludedFiles);

    g_excludedDetails.push_back(exclusion);

  }

}

/**

 * @brief This function excluded frames and logs them

 * Adds excluded frames in the cube files to containers that are used when

 * creating the exclusion file. It also logs the exclusions to the application

 * log window.

 *

 * param listIndex          Location in list for a cube (zero-indexed)

 * param frame              Location in cube for frame (zero-indexed)

 * param reason             Exclusion reason

 *

 * @author 2016-08-17 Victor Silva

 *

 * @internal

 *   @history 2016-08-17 Victor Silva - New app for lroc adapted from makeflat app

 */

void exclude(int listIndex, int frame, string reason){

  if(g_excludedFrames.find(pair<int,int>(listIndex,frame)) == g_excludedFrames.end()) {

    g_excludedFrames.insert(pair < pair < int, int > , string > (pair<int, int>(listIndex, frame), reason));

    PvlObject exclusion("Excluded_Items");

    PvlGroup excludedFiles("Excluded_Frames");

    excludedFiles += PvlKeyword("Frame_from_cube", g_list[listIndex].toString());

    excludedFiles += PvlKeyword("Frame_number", toString(frame));

    excludedFiles += PvlKeyword("Exclusion_reason", QString::fromStdString(reason));

    Application::Log(excludedFiles);

    exclusion.addGroup(excludedFiles);

    g_excludedDetails.push_back(exclusion);

  }

}

BLANKS = "%-6s"    

LENGTH = "%-40s"

include $(ISISROOT)/make/isismake.tststree

APPNAME = lromakeflat

include $(ISISROOT)/make/isismake.tsts

commands:

	ls $(INPUT)/*.cub > $(OUTPUT)/framing.lis

	$(APPNAME) fromlist=$(OUTPUT)/framing.lis \

	  to=$(OUTPUT)/framing \

	  imagetype=framing stdevtol=3 \

	  > /dev/null;

	rm $(OUTPUT)/framing.lis