Statistics.cpp

Go to the documentation of this file.

#include "Statistics.h"

#include <QDataStream>
#include <QDebug>
#include <QString>
#include <QUuid>
#include <QXmlStreamWriter>

#include <float.h>

#include <H5Cpp.h>
#include <hdf5_hl.h>
#include <hdf5.h>

#include "IException.h"
#include "IString.h"
#include "Project.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
#include "XmlStackedHandlerReader.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(Project *project, XmlStackedHandlerReader *xmlReader, QObject *parent) {   // TODO: does xml stuff need project???
//    m_id = NULL;
    SetValidRange();
    Reset(); // initialize
    xmlReader->pushContentHandler(new XmlHandler(this, project));   // TODO: does xml stuff need project???
  }

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

  }


  Statistics::XmlHandler::XmlHandler(Statistics *statistics, Project *project) {   // TODO: does xml stuff need project???
    m_xmlHandlerStatistics = statistics;
    m_xmlHandlerProject = project;   // TODO: does xml stuff need project???
    m_xmlHandlerCharacters = "";
  }


  Statistics::XmlHandler::~XmlHandler() {
    // do not delete this pointer... we don't own it, do we??? passed into StatCumProbDistDynCalc constructor as pointer
    // delete m_xmlHandlerProject;    // TODO: does xml stuff need project???
    m_xmlHandlerProject = NULL;
  }


  bool Statistics::XmlHandler::startElement(const QString &namespaceURI, 
                                                                const QString &localName,
                                                                const QString &qName,
                                                                const QXmlAttributes &atts) {
    m_xmlHandlerCharacters = "";
    if (XmlStackedHandler::startElement(namespaceURI, localName, qName, atts)) {
      // no element attibutes to evaluate
    }
    return true;
  }


  bool Statistics::XmlHandler::characters(const QString &ch) {
    m_xmlHandlerCharacters += ch;
    return XmlStackedHandler::characters(ch);
  }


  bool Statistics::XmlHandler::endElement(const QString &namespaceURI, const QString &localName,
                                     const QString &qName) {
    if (!m_xmlHandlerCharacters.isEmpty()) {
      if (localName == "id") {
//        m_xmlHandlerStatistics->m_id = NULL;
//        m_xmlHandlerStatistics->m_id = new QUuid(m_xmlHandlerCharacters);
      }
      if (localName == "sum") {
        m_xmlHandlerStatistics->m_sum = toDouble(m_xmlHandlerCharacters);
      }
      if (localName == "sumSquares") {
        m_xmlHandlerStatistics->m_sumsum = toDouble(m_xmlHandlerCharacters);
      }
      if (localName == "minimum") {
        m_xmlHandlerStatistics->m_minimum = toDouble(m_xmlHandlerCharacters);
      }
      if (localName == "maximum") {
        m_xmlHandlerStatistics->m_maximum = toDouble(m_xmlHandlerCharacters);
      }
      if (localName == "validMinimum") {
        m_xmlHandlerStatistics->m_validMinimum = toDouble(m_xmlHandlerCharacters);
      }
      if (localName == "validMaximum") {
        m_xmlHandlerStatistics->m_validMaximum = toDouble(m_xmlHandlerCharacters);
      }
      if (localName == "totalPixels") {
        m_xmlHandlerStatistics->m_totalPixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "validPixels") {
        m_xmlHandlerStatistics->m_validPixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "nullPixels") {
        m_xmlHandlerStatistics->m_nullPixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "lisPixels") {
        m_xmlHandlerStatistics->m_lisPixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "lrsPixels") {
        m_xmlHandlerStatistics->m_lrsPixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "hisPixels") {
        m_xmlHandlerStatistics->m_hisPixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "hrsPixels") {
        m_xmlHandlerStatistics->m_hrsPixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "underRangePixels") {
        m_xmlHandlerStatistics->m_underRangePixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "overRangePixels") {
        m_xmlHandlerStatistics->m_overRangePixels = toBigInt(m_xmlHandlerCharacters);
      }
      if (localName == "removedData") {
        m_xmlHandlerStatistics->m_removedData = toBool(m_xmlHandlerCharacters);
      }
      m_xmlHandlerCharacters = "";
    }
    return XmlStackedHandler::endElement(namespaceURI, localName, qName);
  }


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

//??? not working for prog14/15???
// dyld: Symbol not found: __ZN2H58PredType12NATIVE_INT64E
//   Referenced from: /usgs/pkgs/isis3beta2015-12-24/isis/bin/../lib/libisis3.4.12.dylib
//   Expected in: flat namespace
//  in /usgs/pkgs/isis3beta2015-12-24/isis/bin/../lib/libisis3.4.12.dylib

  H5::CompType Statistics::compoundH5DataType() {

    H5::CompType compoundDataType((size_t)124);

    size_t offset = 0;
    compoundDataType.insertMember("Sum", offset, H5::PredType::NATIVE_DOUBLE);

    // offset += sizeof(m_sum);
    offset += sizeof(double);
    compoundDataType.insertMember("SumSquared", offset, H5::PredType::NATIVE_DOUBLE);

    // offset += sizeof(m_sumsum);
    offset += sizeof(double);
    compoundDataType.insertMember("Minimum", offset, H5::PredType::NATIVE_DOUBLE);

    // offset += sizeof(m_minimum);
    offset += sizeof(double);
    compoundDataType.insertMember("Maximum", offset, H5::PredType::NATIVE_DOUBLE);

    // offset += sizeof(m_maximum);
    offset += sizeof(double);
    compoundDataType.insertMember("ValidMinimum", offset, H5::PredType::NATIVE_DOUBLE);

    // offset += sizeof(m_validMinimum);
    offset += sizeof(double);
    compoundDataType.insertMember("ValidMaximum", offset, H5::PredType::NATIVE_DOUBLE);

    // offset += sizeof(m_validMaximum);
    offset += sizeof(double);
    compoundDataType.insertMember("TotalPixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_totalPixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("ValidPixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_validPixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("NullPixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_nullPixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("LRSPixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_lrsPixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("LISPixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_lisPixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("HRSPixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_hrsPixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("HISPixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_hisPixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("UnderRangePixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_underRangePixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("OverRangePixels", offset, H5::PredType::NATIVE_INT64);

    // offset += sizeof(m_overRangePixels);
    offset += sizeof(BigInt);
    compoundDataType.insertMember("RemovedData", offset, H5::PredType::NATIVE_HBOOL);
    // mac osx has problem with "sizeof" commented lines and the native data
    // types too... will consider this later when ready to add serialization.

    return compoundDataType;
  }

} // end namespace isis