Statistics.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
#include "Statistics.h"
 
#include <QDataStream>
#include <QDebug>
#include <QString>
#include <QUuid>
#include <QXmlStreamWriter>
#include <QXmlStreamReader>
 
#include <float.h>
 
#include "IException.h"
#include "IString.h"
#include "Project.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
 
using namespace std;
 
namespace Isis {
  Statistics::Statistics(QObject *parent) : QObject(parent) {
//    m_id = NULL;
//    m_id = new QUuid(QUuid::createUuid());
    SetValidRange();
    Reset(); // initialize
  }
 
  Statistics::Statistics(QXmlStreamReader *xmlReader, QObject *parent) {   // TODO: does xml stuff need project???
//    m_id = NULL;
    SetValidRange();
    Reset(); // initialize
    readStatistics(xmlReader);
  }
 
  void Statistics::readStatistics(QXmlStreamReader *xmlReader) {
    Q_ASSERT(xmlReader->name() == "statistics");
    while (xmlReader->readNextStartElement()) {
      if (xmlReader->qualifiedName() == "sum") {
        try {
          m_sum = toDouble(xmlReader->readElementText());
        }
        catch (IException &e) {
          m_sum = 0.0;
        }
      }
      else if (xmlReader->qualifiedName() == "sumSquares") {
        try {
          m_sumsum = toDouble(xmlReader->readElementText());
        }
        catch (IException &e) {
          m_sumsum = 0.0;
        }
      }
      else if (xmlReader->qualifiedName() == "range") {
        while (xmlReader->readNextStartElement()) {
          if (xmlReader->qualifiedName() == "minimum") {
            try {
              m_minimum = toDouble(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_minimum = DBL_MAX;
            }
          }
          else if (xmlReader->qualifiedName() == "maximum") {
            try {
              m_maximum = toDouble(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_maximum = -DBL_MAX;
            }
          }
          else if (xmlReader->qualifiedName() == "validMinimum") {
            try {
              m_validMinimum = toDouble(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_validMinimum = Isis::ValidMinimum;
            }
          }
          else if (xmlReader->qualifiedName() == "validMaximum") {
            try {
              m_validMaximum = toDouble(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_validMaximum = Isis::ValidMaximum;
            }
          }
          else {
            xmlReader->skipCurrentElement();
          }
        }
      }
      else if (xmlReader->qualifiedName() == "pixelCounts") {
        while (xmlReader->readNextStartElement()) {
          if (xmlReader->qualifiedName() == "totalPixels") {
            try {
              m_totalPixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_totalPixels = 0.0;
            }
          }
          else if (xmlReader->qualifiedName() == "validPixels") {
            try {
              m_validPixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_validPixels = 0.0;
            }
          }
          else if (xmlReader->qualifiedName() == "nullPixels") {
            try {
              m_nullPixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_nullPixels = 0.0;
            }
 
          }
          else if (xmlReader->qualifiedName() == "lisPixels") {
            try {
              m_lisPixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_lisPixels = 0.0;
            }
          }
          else if (xmlReader->qualifiedName() == "lrsPixels") {
            try {
              m_lrsPixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_lrsPixels = 0.0;
            }
          }
          else if (xmlReader->qualifiedName() == "hisPixels") {
            try {
              m_hisPixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_hisPixels = 0.0;
            }
          }
          else if (xmlReader->qualifiedName() == "hrsPixels") {
            try {
              m_hrsPixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_hrsPixels = 0.0;
            }
          }
          else if (xmlReader->qualifiedName() == "underRangePixels") {
            try {
              m_underRangePixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_underRangePixels = 0.0;
            }
          }
          else if (xmlReader->qualifiedName() == "overRangePixels") {
            try {
              m_overRangePixels = toBigInt(xmlReader->readElementText());
            }
            catch (IException &e) {
              m_overRangePixels = 0.0;
            }
          }
          else {
            xmlReader->skipCurrentElement();
          }
        }
      }
      else if (xmlReader->qualifiedName() == "removedData") {
        try {
          m_removedData = toBool(xmlReader->readElementText());
        }
        catch (IException &e) {
          m_removedData = false;
        }
      }
      else {
        xmlReader->skipCurrentElement();
      }
    }
  }
 
  Statistics::Statistics(const PvlGroup &inStats, QObject *parent) {
    SetValidRange();
    Reset();
    fromPvl(inStats);
  }
 
 
  Statistics::Statistics(const Statistics &other)
    : m_sum(other.m_sum),
      m_sumsum(other.m_sumsum),
      m_minimum(other.m_minimum),
      m_maximum(other.m_maximum),
      m_validMinimum(other.m_validMinimum),
      m_validMaximum(other.m_validMaximum),
      m_totalPixels(other.m_totalPixels),
      m_validPixels(other.m_validPixels),
      m_nullPixels(other.m_nullPixels),
      m_lrsPixels(other.m_lrsPixels),
      m_lisPixels(other.m_lisPixels),
      m_hrsPixels(other.m_hrsPixels),
      m_hisPixels(other.m_hisPixels),
      m_underRangePixels(other.m_underRangePixels),
      m_overRangePixels(other.m_overRangePixels),
      m_removedData(other.m_removedData) {
  }
   // : m_id(new QUuid(other.m_id->toString())),
 
 
  Statistics::~Statistics() {
//    delete m_id;
//    m_id = NULL;
  }
 
 
  Statistics &Statistics::operator=(const Statistics &other) {
 
    if (&other != this) {
//      delete m_id;
//      m_id = NULL;
//      m_id = new QUuid(other.m_id->toString());
 
      m_sum = other.m_sum;
      m_sumsum = other.m_sumsum;
      m_minimum = other.m_minimum;
      m_maximum = other.m_maximum;
      m_validMinimum = other.m_validMinimum;
      m_validMaximum = other.m_validMaximum;
      m_totalPixels = other.m_totalPixels;
      m_validPixels = other.m_validPixels;
      m_nullPixels = other.m_nullPixels;
      m_lrsPixels = other.m_lrsPixels;
      m_lisPixels = other.m_lisPixels;
      m_hrsPixels = other.m_hrsPixels;
      m_hisPixels = other.m_hisPixels;
      m_underRangePixels = other.m_underRangePixels;
      m_overRangePixels = other.m_overRangePixels;
      m_removedData = other.m_removedData;
    }
    return *this;
 
  }
 
 
  void Statistics::Reset() {
    m_sum = 0.0;
    m_sumsum = 0.0;
    m_minimum = DBL_MAX;
    m_maximum = -DBL_MAX;
    m_totalPixels = 0;
    m_validPixels = 00;
    m_nullPixels = 0;
    m_lisPixels = 0;
    m_lrsPixels = 0;
    m_hrsPixels = 0;
    m_hisPixels = 0;
    m_overRangePixels = 0;
    m_underRangePixels = 0;
    m_removedData = false;
  }
 
 
  void Statistics::AddData(const double *data, const unsigned int count) {
    for(unsigned int i = 0; i < count; i++) {
      double value = data[i];
      AddData(value);
    }
  }
 
 
  void Statistics::AddData(const double data) {
    m_totalPixels++;
 
    if (Isis::IsNullPixel(data)) {
      m_nullPixels++;
    }
    else if (Isis::IsHisPixel(data)) {
      m_hisPixels++;
    }
    else if (Isis::IsHrsPixel(data)) {
      m_hrsPixels++;
    }
    else if (Isis::IsLisPixel(data)) {
      m_lisPixels++;
    }
    else if (Isis::IsLrsPixel(data)) {
      m_lrsPixels++;
    }
    else if (AboveRange(data)) {
      m_overRangePixels++;
    }
    else if (BelowRange(data)) {
      m_underRangePixels++;
    }
    else { // if (Isis::IsValidPixel(data) && InRange(data)) {
      m_sum += data;
      m_sumsum += data * data;
      if (data < m_minimum) m_minimum = data;
      if (data > m_maximum) m_maximum = data;
      m_validPixels++;
    }
 
  }
 
 
  void Statistics::RemoveData(const double *data, const unsigned int count) {
 
    for(unsigned int i = 0; i < count; i++) {
      double value = data[i];
      RemoveData(value);
    }
 
  }
 
 
  void Statistics::RemoveData(const double data) {
    m_removedData = true;
    m_totalPixels--;
 
    if (Isis::IsNullPixel(data)) {
      m_nullPixels--;
    }
    else if (Isis::IsHisPixel(data)) {
      m_hisPixels--;
    }
    else if (Isis::IsHrsPixel(data)) {
      m_hrsPixels--;
    }
    else if (Isis::IsLisPixel(data)) {
      m_lisPixels--;
    }
    else if (Isis::IsLrsPixel(data)) {
      m_lrsPixels--;
    }
    else if (AboveRange(data)) {
      m_overRangePixels--;
    }
    else if (BelowRange(data)) {
      m_underRangePixels--;
    }
    else { // if (IsValidPixel(data) && InRange(data)) {
      m_sum -= data;
      m_sumsum -= data * data;
      m_validPixels--;
    }
 
    if (m_totalPixels < 0) {
      QString msg = "You are removing non-existant data in [Statistics::RemoveData]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    // what happens to saved off min/max???
  }
 
 
  void Statistics::SetValidRange(const double minimum, const double maximum) {
    m_validMinimum = minimum;
    m_validMaximum = maximum;
 
    if (m_validMaximum < m_validMinimum) {
      // get the min and max DN values in the chosen range
      QString msg = "Invalid Range: Minimum [" + toString(minimum) 
                    + "] must be less than the Maximum [" + toString(maximum) + "].";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
    //??? throw exception if data has already been added???
  }
 
 
  double Statistics::ValidMinimum() const {
    return m_validMinimum;
  }
 
 
  double Statistics::ValidMaximum() const {
    return m_validMaximum;
  }
 
 
  bool Statistics::InRange(const double value) {
    return (!BelowRange(value) && !AboveRange(value));
  }
 
 
  bool Statistics::AboveRange(const double value) {
    return (value > m_validMaximum);
  }
 
 
  bool Statistics::BelowRange(const double value) {
    return (value < m_validMinimum);
  }
 
 
  double Statistics::Average() const {
    if (m_validPixels < 1) return Isis::NULL8;
    return m_sum / m_validPixels;
  }
 
 
  double Statistics::StandardDeviation() const {
    if (m_validPixels <= 1) return Isis::NULL8;
    return sqrt(Variance());
  }
 
 
  double Statistics::Variance() const {
    if (m_validPixels <= 1) return Isis::NULL8;
    double temp = m_validPixels * m_sumsum - m_sum * m_sum;
    if (temp < 0.0) temp = 0.0;  // This should happen unless roundoff occurs
    return temp / ((m_validPixels - 1.0) * m_validPixels);
  }
 
 
  double Statistics::Sum() const {
    return m_sum;
  }
 
 
  double Statistics::SumSquare() const {
    return m_sumsum;
  }
 
 
  double Statistics::Rms() const {
    if (m_validPixels < 1) return Isis::NULL8;
    double temp = m_sumsum / m_validPixels;
    if (temp < 0.0) temp = 0.0;
    return sqrt(temp);
  }
 
 
  double Statistics::Minimum() const {
    if (m_removedData) {
      QString msg = "Minimum is invalid since you removed data";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if (m_validPixels < 1) return Isis::NULL8;
    return m_minimum;
  }
 
 
  double Statistics::Maximum() const {
    if (m_removedData) {
      QString msg = "Maximum is invalid since you removed data";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if (m_validPixels < 1) return Isis::NULL8;
    return m_maximum;
  }
 
 
  BigInt Statistics::TotalPixels() const {
    return m_totalPixels;
  }
 
 
  BigInt Statistics::ValidPixels() const {
    return m_validPixels;
  }
 
 
  BigInt Statistics::OverRangePixels() const {
    return m_overRangePixels;
  }
 
 
  BigInt Statistics::UnderRangePixels() const {
    return m_underRangePixels;
  }
 
 
  BigInt Statistics::NullPixels() const {
    return m_nullPixels;
  }
 
 
  BigInt Statistics::LisPixels() const {
    return m_lisPixels;
  }
 
 
  BigInt Statistics::LrsPixels() const {
    return m_lrsPixels;
  }
 
 
  BigInt Statistics::HisPixels() const {
    return m_hisPixels;
  }
 
 
  BigInt Statistics::HrsPixels() const {
    return m_hrsPixels;
  }
 
 
  BigInt Statistics::OutOfRangePixels() const {
    return m_overRangePixels + m_underRangePixels;
  }
 
 
  bool Statistics::RemovedData() const {
    return m_removedData;
  }
 
 
  double Statistics::ChebyshevMinimum(const double percent) const {
    if ((percent <= 0.0) || (percent >= 100.0)) {
      QString msg = "Invalid value for percent";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if (m_validPixels < 1) return Isis::NULL8;
    double k = sqrt(1.0 / (1.0 - percent / 100.0));
    return Average() - k * StandardDeviation();
  }
 
 
  double Statistics::ChebyshevMaximum(const double percent) const {
    if ((percent <= 0.0) || (percent >= 100.0)) {
      QString msg = "Invalid value for percent";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }
 
    if (m_validPixels < 1) return Isis::NULL8;
    double k = sqrt(1.0 / (1.0 - percent / 100.0));
    return Average() + k * StandardDeviation();
  }
 
 
  double Statistics::BestMinimum(const double percent) const {
    if (m_validPixels < 1) return Isis::NULL8;
    // ChebyshevMinimum can return erroneous values when the data is all a
    // single value
    // In this case, we just want to return the minimum
    if (Minimum() == Maximum()) return Minimum();
    double min = ChebyshevMinimum(percent);
    if (Minimum() > min) min = Minimum();
    return min;
  }
 
 
  double Statistics::BestMaximum(const double percent) const {
    if (m_validPixels < 1) return Isis::NULL8;
    // ChebyshevMaximum can return erroneous values when the data is all a
    // single value
    // In this case, we just want to return the maximum
    if (Minimum() == Maximum()) return Maximum();
    double max = ChebyshevMaximum(percent);
    if (Maximum() < max) max = Maximum();
    return max;
  }
 
 
  double Statistics::ZScore(const double value) const {
    if (StandardDeviation() == 0) {
      if (value == Maximum()) return 0;
      else {
        QString msg = "Undefined Z-score. Standard deviation is zero and the input value[" 
                      + toString(value) + "] is out of range [" + toString(Maximum()) + "].";
        throw IException(IException::Programmer, msg, _FILEINFO_);
      }
    }
    return (value - Average()) / StandardDeviation();
  }
 
 
  void Statistics::fromPvl(const PvlGroup &inStats) {
    Reset();
    m_sum = inStats["Sum"];
    m_sumsum = inStats["SumSquare"];
    m_minimum = inStats["Minimum"];
    m_maximum = inStats["Maximum"];
    m_validMinimum = inStats["ValidMinimum"];
    m_validMaximum = inStats["ValidMaximum"];
    m_totalPixels = inStats["TotalPixels"];
    m_validPixels = inStats["ValidPixels"];
    m_nullPixels = inStats["NullPixels"];
    m_lrsPixels = inStats["LrsPixels"];
    m_lisPixels = inStats["LisPixels"];
    m_hrsPixels = inStats["HrsPixels"];
    m_hisPixels = inStats["HisPixels"];
    m_underRangePixels = inStats["UnderValidMinimumPixels"];
    m_overRangePixels = inStats["OverValidMaximumPixels"];
    m_removedData = false; //< Is this the correct state?
  }
 
 
  PvlGroup Statistics::toPvl(QString name) const {
    if (name.isEmpty()) {
      name = "Statistics";
    }
    // Construct a label with the results
    PvlGroup results(name);  
    results += PvlKeyword("Sum", toString(Sum()));
    results += PvlKeyword("SumSquare", toString(SumSquare()));
    results += PvlKeyword("Minimum", toString(Minimum()));
    results += PvlKeyword("Maximum", toString(Maximum()));
    results += PvlKeyword("ValidMinimum", toString(ValidMinimum()));
    results += PvlKeyword("ValidMaximum", toString(ValidMaximum()));
    if (ValidPixels() != 0) {
      results += PvlKeyword("Average", toString(Average()));
      results += PvlKeyword("StandardDeviation", toString(StandardDeviation()));
      results += PvlKeyword("Variance", toString(Variance()));
    }
    results += PvlKeyword("TotalPixels", toString(TotalPixels()));
    results += PvlKeyword("ValidPixels", toString(ValidPixels()));
    results += PvlKeyword("OverValidMaximumPixels", toString(OverRangePixels()));
    results += PvlKeyword("UnderValidMinimumPixels", toString(UnderRangePixels()));
    results += PvlKeyword("NullPixels", toString(NullPixels()));
    results += PvlKeyword("LisPixels", toString(LisPixels()));
    results += PvlKeyword("LrsPixels", toString(LrsPixels()));
    results += PvlKeyword("HisPixels", toString(HisPixels()));
    results += PvlKeyword("HrsPixels", toString(HrsPixels()));
 
    return results;
  }
 
 
  void Statistics::save(QXmlStreamWriter &stream, const Project *project) const {   // TODO: does xml stuff need project???
 
    stream.writeStartElement("statistics");
//    stream.writeTextElement("id", m_id->toString());
 
    stream.writeTextElement("sum", toString(m_sum));
    stream.writeTextElement("sumSquares", toString(m_sumsum));
 
    stream.writeStartElement("range");
    stream.writeTextElement("minimum", toString(m_minimum));
    stream.writeTextElement("maximum", toString(m_maximum));
    stream.writeTextElement("validMinimum", toString(m_validMinimum));
    stream.writeTextElement("validMaximum", toString(m_validMaximum));
    stream.writeEndElement(); // end range
    
    stream.writeStartElement("pixelCounts");
    stream.writeTextElement("totalPixels", toString(m_totalPixels));
    stream.writeTextElement("validPixels", toString(m_validPixels));
    stream.writeTextElement("nullPixels", toString(m_nullPixels));
    stream.writeTextElement("lisPixels", toString(m_lisPixels));
    stream.writeTextElement("lrsPixels", toString(m_lrsPixels));
    stream.writeTextElement("hisPixels", toString(m_hisPixels));
    stream.writeTextElement("hrsPixels", toString(m_hrsPixels));
    stream.writeTextElement("underRangePixels", toString(m_underRangePixels));
    stream.writeTextElement("overRangePixels", toString(m_overRangePixels));
    stream.writeEndElement(); // end pixelCounts
    
    stream.writeTextElement("removedData", toString(m_removedData));
    stream.writeEndElement(); // end statistics
 
  }
 
 
  QDataStream &Statistics::write(QDataStream &stream) const {
    stream.setByteOrder(QDataStream::LittleEndian);
    stream << m_sum
           << m_sumsum
           << m_minimum
           << m_maximum
           << m_validMinimum
           << m_validMaximum
           << (qint64)m_totalPixels
           << (qint64)m_validPixels
           << (qint64)m_nullPixels
           << (qint64)m_lrsPixels
           << (qint64)m_lisPixels
           << (qint64)m_hrsPixels
           << (qint64)m_hisPixels
           << (qint64)m_underRangePixels
           << (qint64)m_overRangePixels
           << (qint32)m_removedData;
    return stream;
//    stream << m_id->toString()
  }
 
 
  QDataStream &Statistics::read(QDataStream &stream) {
 
//    QString id;
    qint64 totalPixels, validPixels, nullPixels, lrsPixels, lisPixels,
           hrsPixels, hisPixels, underRangePixels, overRangePixels;
    qint32 removedData;
 
    stream.setByteOrder(QDataStream::LittleEndian);
//    stream >> id
    stream >> m_sum
           >> m_sumsum
           >> m_minimum
           >> m_maximum
           >> m_validMinimum
           >> m_validMaximum
           >> totalPixels
           >> validPixels
           >> nullPixels
           >> lrsPixels
           >> lisPixels
           >> hrsPixels
           >> hisPixels
           >> underRangePixels
           >> overRangePixels
           >> removedData;
 
//    delete m_id;
//    m_id = NULL;
//    m_id = new QUuid(id);
 
    m_totalPixels      = (BigInt)totalPixels;
    m_validPixels      = (BigInt)validPixels;
    m_nullPixels       = (BigInt)nullPixels;
    m_lrsPixels        = (BigInt)lrsPixels;
    m_lisPixels        = (BigInt)lisPixels;
    m_hrsPixels        = (BigInt)hrsPixels;
    m_hisPixels        = (BigInt)hisPixels;
    m_underRangePixels = (BigInt)underRangePixels;
    m_overRangePixels  = (BigInt)overRangePixels;
    m_removedData      = (bool)removedData;
    
    return stream;
  }
 
 
  QDataStream &operator<<(QDataStream &stream, const Statistics &statistics) {
    return statistics.write(stream);
  }
 
 
  QDataStream &operator>>(QDataStream &stream, Statistics &statistics) {
    return statistics.read(stream);
  }
 
} // end namespace isis