Histogram.cpp

Go to the documentation of this file.

#include "Histogram.h"

#include <string>
#include <iostream>
#include <stdio.h>

#include "Brick.h"
#include "ControlMeasure.h"
#include "ControlNet.h"
#include "ControlPoint.h"
#include "Message.h"
#include "LineManager.h"
#include <math.h>

using namespace std;
namespace Isis {

  Histogram::Histogram(double minimum, double maximum, int nbins) {
    SetValidRange(minimum, maximum);
    SetBinRange(minimum, maximum);
    SetBins(nbins);
  }


  Histogram::Histogram(Cube &cube, int statsBand, Progress *progress,
      double startSample, double startLine,
      double endSample, double endLine,
      int bins, bool addCubeData) {
    InitializeFromCube(cube, statsBand, progress, bins, startSample, startLine,
                       endSample, endLine);

    if (addCubeData) {
      Brick cubeDataBrick((int)(endSample - startSample + 1),
                          1, 1, cube.pixelType());

      // if band == 0, then we're gathering data for all bands.
      int startBand = statsBand;
      int endBand = statsBand;

      if (statsBand == 0) {
        startBand = 1;
        endBand = cube.bandCount();
      }

      if (progress != NULL) {
        progress->SetText("Gathering histogram");
        progress->SetMaximumSteps(
          (int)(endLine - startLine + 1) * (int)(endBand - startBand + 1));
        progress->CheckStatus();
      }

      for (int band = startBand; band <= endBand; band++) {
        for (int line = (int)startLine; line <= endLine; line++) {
          cubeDataBrick.SetBasePosition(qRound(startSample), line, band);
          cube.read(cubeDataBrick);
          AddData(cubeDataBrick.DoubleBuffer(), cubeDataBrick.size());
          if (progress != NULL) {
            progress->CheckStatus();
          }
        }
      }
    }
  }


  Histogram::Histogram(ControlNet &net, double(ControlMeasure::*statFunc)() const, int bins) {
    
    //check to make sure we have a reasonable number of bins
    if (bins < 1) {
      string msg = "The number of Histogram Bins must be greater than 0";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    SetBins(bins);

    //set the ranges
    rangesFromNet(net,statFunc);

    //add all the data to the now setup histogram
    addMeasureDataFromNet(net,statFunc);
  }


  Histogram::Histogram(ControlNet &net, double(ControlMeasure::*statFunc)() const, double binWidth) {
    
    //check to make sure we have a reasonable number of bins
    if (binWidth <= 0 ) {
      string msg = "The width of Histogram Bins must be greater than 0";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

   //get the range of the data
   rangesFromNet(net,statFunc);

   
   //stretch the domain so that it is an even multiple of binWidth
     //for some reason Histogram makes the end points of the bin range be at the center of
     //bins.  Thus the +/-0.5 forces it to point the bin range at the ends of the bins.
   SetBinRange(binWidth*(floor(this->ValidMinimum()/binWidth)+0.5),
               binWidth*(ceil(this->ValidMaximum()/binWidth)-0.5));
   SetValidRange(binWidth*floor(this->ValidMinimum()/binWidth),
                 binWidth*ceil(this->ValidMaximum()/binWidth));

   //from the domain of the data and the requested bin width calculate the number of bins
   double domain = this->ValidMaximum() - this->ValidMinimum();
   int nBins = int ( ceil(domain/binWidth) );
   SetBins(nBins);

   //add all the data to the now setup histogram
   addMeasureDataFromNet(net,statFunc);
  }


  void Histogram::addMeasureDataFromNet(ControlNet &net, double(ControlMeasure::*statFunc)() const) {
    //get the number of object points
    int nObjPts =  net.GetNumPoints();
    for (int i=0;i<nObjPts;i++) { //for each Object point
      const ControlPoint *point = net.GetPoint(i);
      if (point->IsIgnored()) continue;  //if the point is ignored then continue

      //get the number of measures
      int nObs = point->GetNumMeasures();
      for (int j=0;j<nObs;j++) {  //for every measure
        const ControlMeasure *measure = point->GetMeasure(j);
        if (measure->IsIgnored())  continue;
        this->AddData((measure->*statFunc)());
      }
    }
  }


  void Histogram::rangesFromNet(ControlNet &net, double(ControlMeasure::*statFunc)() const) {
    double min,max,temp;
    min =  DBL_MAX;
    max = -DBL_MAX;

    //get the number of object points
    int nObjPts =  net.GetNumPoints();
    for (int i=0;i<nObjPts;i++) { //for each Object point
      const ControlPoint *point = net.GetPoint(i);
      if (point->IsIgnored()) continue;  //if the point is ignored then continue

      //get the number of measures
      int nObs = point->GetNumMeasures();
      for (int j=0;j<nObs;j++) {  //for every measure
        const ControlMeasure *measure = point->GetMeasure(j);
        if (measure->IsIgnored())  continue;

        temp = (measure->*statFunc)();  //get the data using the passed ControlMeasure acessor function pointer
        if (temp > max) max = temp;
        if (temp < min) min = temp;
      }
    }

    //if DBL_MAX's weren't changed there's a problem
    if (max <= min) {
      string msg = "The net file appears to have 1 or fewer measures with residual data, "
                   "thus no histogram for this net file can be created;";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    //set up the histogram ranges
    SetValidRange(min, max);
    SetBinRange(min, max);
  }


  void Histogram::InitializeFromCube(Cube &cube, int statsBand,
      Progress *progress, int nbins, double startSample, double startLine,
      double endSample, double endLine) {
    // Make sure band is valid, 0 is valid (means all bands)
    if ((statsBand < 0) || (statsBand > cube.bandCount())) {
      string msg = "Cannot gather histogram for band [" + IString(statsBand) +
          "]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // We need to find the min/max DN value for our histogram bins to be the
    //   correct size.
    double minDnValue = Null;
    double maxDnValue = Null;

    if (cube.pixelType() == UnsignedByte) {
      // If we can discretely store every data point, then we can use the
      //   possible extent of the data range as our min/max dn values.
      if (nbins == 0) {
        minDnValue = 0.0 * cube.multiplier() + cube.base();
        maxDnValue = 255.0 * cube.multiplier() + cube.base();
        nbins = 256;
      }
    }
    else if (cube.pixelType() == SignedWord) {
      if (nbins == 0) {
        minDnValue = -32768.0 * cube.multiplier() + cube.base();
        maxDnValue = 32767.0 * cube.multiplier() + cube.base();
        nbins = 65536;
      }
    }
    else if (cube.pixelType() == Real) {
      // We can't account for all the possibilities of a double inside of any
      //   data range, so don't guess min/max DN values.
      if (nbins == 0) {
        nbins = 65536;
      }
    }
    else {
      IString msg = "Unsupported pixel type";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if (startSample == Null)
      startSample = 1.0;

    if (endSample == Null)
      endSample = cube.sampleCount();

    if (startLine == Null)
      startLine = 1.0;

    if (endLine == Null)
      endLine = cube.lineCount();

    // If we still need our min/max DN values, find them.
    if (minDnValue == Null || maxDnValue == Null) {

      Brick cubeDataBrick((int)(endSample - startSample + 1),
                          1, 1, cube.pixelType());
      Statistics stats;

      // if band == 0, then we're gathering stats for all bands. I'm really
      //   not sure that this is correct, a good idea or called from anywhere...
      //   but I don't have time to track down the use case.
      int startBand = statsBand;
      int endBand = statsBand;

      if (statsBand == 0) {
        startBand = 1;
        endBand = cube.bandCount();
      }

      if (progress != NULL) {
        progress->SetText("Computing min/max for histogram");
        progress->SetMaximumSteps(
          (int)(endLine - startLine + 1) * (int)(endBand - startBand + 1));
        progress->CheckStatus();
      }

      for (int band = startBand; band <= endBand; band++) {
        for (int line = (int)startLine; line <= endLine; line++) {
          cubeDataBrick.SetBasePosition(qRound(startSample), line, band);
          cube.read(cubeDataBrick);
          stats.AddData(cubeDataBrick.DoubleBuffer(), cubeDataBrick.size());
          if (progress != NULL) {
            progress->CheckStatus();
          }
        }
      }

      if (stats.ValidPixels() == 0) {
        minDnValue = 0.0;
        maxDnValue = 1.0;
      }
      else {
        minDnValue = stats.Minimum();
        maxDnValue = stats.Maximum();
      }
    }

    // Set the bins and range
    SetBinRange(minDnValue, maxDnValue);
    SetBins(nbins);
  }


  Histogram::~Histogram() {
  }

  void Histogram::SetBinRange(double binStart, double binEnd) {
    if (binEnd < binStart) {
      string msg = "The binning range start [" + IString(binStart) +
                   " must be less than the end [" + IString(binEnd) + ".";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    p_binRangeStart = binStart;
    p_binRangeEnd = binEnd;
  }


  void Histogram::Reset() {
    Statistics::Reset();
    for (int i = 0; i < (int)p_bins.size(); i++) {
      p_bins[i] = 0;
    }
  }


  void Histogram::SetBins(const int nbins) {
    p_bins.resize(nbins);
    Histogram::Reset();
  }

  void Histogram::AddData(const double *data,
                          const unsigned int count) {
    Statistics::AddData(data, count);

    int nbins = p_bins.size();
    int index;
    for (unsigned int i = 0; i < count; i++) {
      if (IsValidPixel(data[i]) && InRange(data[i])) {
        if (BinRangeStart() == BinRangeEnd()) {
          index = 0;
        }
        else {
          index = (int) floor((double)(nbins - 1) / (BinRangeEnd() - BinRangeStart()) *
                              (data[i] - BinRangeStart()) + 0.5);
        }
        if (index < 0) index = 0;
        if (index >= nbins) index = nbins - 1;
        p_bins[index] += 1;
      }
    }
  }


  void Histogram::AddData(const double data) {
    Statistics::AddData(data);

    int nbins = p_bins.size();
    int index;
    if (IsValidPixel(data) && InRange(data)) {
      if (BinRangeStart() == BinRangeEnd()) {
        index = 0;
      }
      else {
        index = (int) floor((double)(nbins - 1) / (BinRangeEnd() - BinRangeStart()) *
                            (data - BinRangeStart()) + 0.5);
      }
      if (index < 0) index = 0;
      if (index >= nbins) index = nbins - 1;
      p_bins[index] += 1;
    }
  }


  void Histogram::RemoveData(const double *data,
                             const unsigned int count) {
    Statistics::RemoveData(data, count);

    int nbins = p_bins.size();
    int index;
    for (unsigned int i = 0; i < count; i++) {
      if (IsValidPixel(data[i])) {
        if (BinRangeStart() == BinRangeEnd()) {
          index = 0;
        }
        else {
          index = (int) floor((double)(nbins - 1) / (BinRangeEnd() - BinRangeStart()) *
                              (data[i] - BinRangeStart()) + 0.5);
        }
        if (index < 0) index = 0;
        if (index >= nbins) index = nbins - 1;
        p_bins[index] -= 1;
      }
    }
  }

  double Histogram::Median() const {
    return Percent(50.0);
  }

  double Histogram::Mode() const {
    int mode = 0;
    for (int i = 0; i < (int)p_bins.size(); i++) {
      if (p_bins[i] > p_bins[mode]) mode = i;
    }

    if (p_bins[mode] < 1) return NULL8;

    return BinMiddle(mode);
  }


  double Histogram::Percent(double percent) const {
    if ((percent < 0.0) || (percent > 100.0)) {
      string m = "Argument percent outside of the range 0 to 100 in";
      m += " [Histogram::Percent]";
      throw IException(IException::Programmer, m, _FILEINFO_);
    }

    if (ValidPixels() < 1) return NULL8;

    BigInt currentPixels = 0;
    double currentPercent;

    for (int i = 0; i < (int)p_bins.size(); i++) {
      currentPixels += p_bins[i];
      currentPercent = (double) currentPixels / (double) ValidPixels() * 100.0;
      if (currentPercent >= percent) {
        return BinMiddle(i);
      }
    }

    return BinMiddle((int)p_bins.size() - 1);
  }


  double Histogram::Skew() const {
    if (ValidPixels() < 1) return NULL8;
    double sdev = StandardDeviation();
    if (sdev == 0.0) return 0.0;
    return 3.0 * (Average() - Median()) / sdev;
  }


  BigInt Histogram::BinCount(const int index) const {
    if ((index < 0) || (index >= (int)p_bins.size())) {
      QString message = Message::ArraySubscriptNotInRange(index);
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    return p_bins[index];
  }


  void Histogram::BinRange(const int index,
                           double &low, double &high) const {
    if ((index < 0) || (index >= (int)p_bins.size())) {
      QString message = Message::ArraySubscriptNotInRange(index);
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    double binSize = (BinRangeEnd() - BinRangeStart()) / (double)(p_bins.size() - 1);
    low = BinRangeStart() - binSize / 2.0 + binSize * (double) index;
    high = low + binSize;
  }


  double Histogram::BinMiddle(const int index) const {
    if ((index < 0) || (index >= (int)p_bins.size())) {
      QString message = Message::ArraySubscriptNotInRange(index);
      throw IException(IException::Programmer, message, _FILEINFO_);
    }

    double low, high;
    BinRange(index, low, high);
    return (low + high) / 2.0;
  }


  double Histogram::BinSize() const {
    double low, high;
    BinRange(0, low, high);
    return high - low;
  }


  int Histogram::Bins() const {
    return (int)p_bins.size();
  }


  BigInt Histogram::MaxBinCount() const {

    BigInt maxBinCount = 0;
    for (int i = 0; i < (int)p_bins.size(); i++) {
      if (p_bins[i] > maxBinCount) maxBinCount = p_bins[i];
    }

    return maxBinCount;
  }
}