Histogram.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "Histogram.h"
 
#include "ControlMeasure.h"
#include "ControlNet.h"
#include "ControlPoint.h"
#include "LineManager.h"
#include "Message.h"
#include "SpecialPixel.h"
 
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <string>
 
using namespace std;
 
namespace Isis {
 
  Histogram::Histogram(double minimum, double maximum, int nbins) {
 
    SetValidRange(minimum, maximum);
    SetBins(nbins);
  }
 
 
  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);
 
   //Keep an eye on this to see if it breaks anything.  Also, I need to create
   //a dataset to give this constructor for the unit test.
 
   //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) ) - 1;
   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 ( !IsSpecial(temp) && temp > max ) max = temp;
        if ( !IsSpecial(temp) && 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);
  }
 
 
  Histogram::~Histogram() {
  }
 
  //2015-08-24,  Tyler Wilson:  Added Statistics::SetValidRange call to SetValidRange
  //So the two functions do not have to be called together when setting
  //up a histogram
 
  void Histogram::SetValidRange(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_);
    }
 
    //(tjw): Flush the data buffers.  Since we are setting
    //the statistical range for the data, any data loaded
    //before this call is useless.
    Reset();
    Isis::Statistics::SetValidRange(binStart,binEnd);
    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 / (BinRangeEnd() - BinRangeStart())) *
                              (data[i] - BinRangeStart()) );
        }
        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 / (BinRangeEnd() - BinRangeStart())) *
                              (data - BinRangeStart()) );
      }
      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 / (BinRangeEnd() - BinRangeStart())) *
                              (data[i] - BinRangeStart()) );
        }
        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();
    low = BinRangeStart() + 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;
  }
}