OverlapStatistics.cpp

Go to the documentation of this file.

#include "OverlapStatistics.h"

#include <cfloat>
#include <iomanip>

#include "Brick.h"
#include "Cube.h"
#include "FileName.h"
#include "IException.h"
#include "MultivariateStatistics.h"
#include "Progress.h"
#include "Projection.h"
#include "ProjectionFactory.h"
#include "PvlGroup.h"
#include "PvlKeyword.h"
#include "PvlObject.h"

using namespace std;
namespace Isis {

  OverlapStatistics::OverlapStatistics(const PvlObject &inStats) {
    init();
    fromPvl(inStats);
  }

  OverlapStatistics::OverlapStatistics(Isis::Cube &x, Isis::Cube &y,
                                       QString progressMsg, double sampPercent) {

    init();

    // Test to ensure sampling percent in bound
    if (sampPercent <= 0.0 || sampPercent > 100.0) {
      string msg = "The sampling percent must be a decimal (0.0, 100.0]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    p_sampPercent = sampPercent;

    // Extract filenames and band number from cubes
    p_xFile = x.fileName();
    p_yFile = y.fileName();

    // Make sure number of bands match
    if (x.bandCount() != y.bandCount()) {
      QString msg = "Number of bands do not match between cubes [" +
                   p_xFile.name() + "] and [" + p_yFile.name() + "]";
      throw IException(IException::User, msg, _FILEINFO_);
    }
    p_bands = x.bandCount();
    p_stats.resize(p_bands);

    //Create projection from each cube
    Projection *projX = x.projection();
    Projection *projY = y.projection();

    // Test to make sure projection parameters match
    if (*projX != *projY) {
      QString msg = "Mapping groups do not match between cubes [" +
                   p_xFile.name() + "] and [" + p_yFile.name() + "]";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    // Figure out the x/y range for both images to find the overlap
    double Xmin1 = projX->ToProjectionX(0.5);
    double Ymax1 = projX->ToProjectionY(0.5);
    double Xmax1 = projX->ToProjectionX(x.sampleCount() + 0.5);
    double Ymin1 = projX->ToProjectionY(x.lineCount() + 0.5);

    double Xmin2 = projY->ToProjectionX(0.5);
    double Ymax2 = projY->ToProjectionY(0.5);
    double Xmax2 = projY->ToProjectionX(y.sampleCount() + 0.5);
    double Ymin2 = projY->ToProjectionY(y.lineCount() + 0.5);

    // Find overlap
    if ((Xmin1 < Xmax2) && (Xmax1 > Xmin2) && (Ymin1 < Ymax2) && (Ymax1 > Ymin2)) {
      double minX = Xmin1 > Xmin2 ? Xmin1 : Xmin2;
      double minY = Ymin1 > Ymin2 ? Ymin1 : Ymin2;
      double maxX = Xmax1 < Xmax2 ? Xmax1 : Xmax2;
      double maxY = Ymax1 < Ymax2 ? Ymax1 : Ymax2;

      // Find Sample range of the overlap
      p_minSampX = (int)(projX->ToWorldX(minX) + 0.5);
      p_maxSampX = (int)(projX->ToWorldX(maxX) + 0.5);
      p_minSampY = (int)(projY->ToWorldX(minX) + 0.5);
      p_maxSampY = (int)(projY->ToWorldX(maxX) + 0.5);
      p_sampRange = p_maxSampX - p_minSampX + 1;

      // Test to see if there was only sub-pixel overlap
      if (p_sampRange <= 0) return;

      // Find Line range of overlap
      p_minLineX = (int)(projX->ToWorldY(maxY) + 0.5);
      p_maxLineX = (int)(projX->ToWorldY(minY) + 0.5);
      p_minLineY = (int)(projY->ToWorldY(maxY) + 0.5);
      p_maxLineY = (int)(projY->ToWorldY(minY) + 0.5);
      p_lineRange = p_maxLineX - p_minLineX + 1;

      // Print percent processed
      Progress progress;
      progress.SetText(progressMsg);

      int linc = (int)(100.0 / sampPercent + 0.5); // Calculate our line increment

      // Define the maximum number of steps to be our line range divided by the
      // line increment, but if they do not divide evenly, then because of
      // rounding, we need to do an additional step for each band
      int maxSteps = (int)(p_lineRange / linc + 0.5);

      if (p_lineRange % linc != 0) maxSteps += 1;
      maxSteps *= p_bands;


      progress.SetMaximumSteps(maxSteps);
      progress.CheckStatus();

      // Collect and store off the overlap statistics
      for (int band = 1; band <= p_bands; band++) {
        Brick b1(p_sampRange, 1, 1, x.pixelType());
        Brick b2(p_sampRange, 1, 1, y.pixelType());

        int i = 0;
        while(i < p_lineRange) {
          b1.SetBasePosition(p_minSampX, (i + p_minLineX), band);
          b2.SetBasePosition(p_minSampY, (i + p_minLineY), band);
          x.read(b1);
          y.read(b2);
          p_stats[band-1].AddData(b1.DoubleBuffer(), b2.DoubleBuffer(), p_sampRange);

          // Make sure we consider the last line
          if (i + linc > p_lineRange - 1 && i != p_lineRange - 1) {
            i = p_lineRange - 1;
            progress.AddSteps(1);
          }
          else i += linc; // Increment the current line by our incrementer

          progress.CheckStatus();
        }
      }
    }
  }


  bool OverlapStatistics::HasOverlap() const {
    for (int b = 0; b < p_bands; b++) {
      if (p_stats[b].ValidPixels() > 0) return true;
    }
    return false;
  }


  PvlObject OverlapStatistics::toPvl(QString name) const {
    try {
      // Empty strings are set to default value as well
      if (name.isEmpty()) {
        name = "OverlapStatistics";
      }

      // Add keywords for OverlapStatistics data
      PvlObject o(name);
      o += PvlKeyword("File1", FileNameX().name());
      o += PvlKeyword("File2", FileNameY().name());
      o += PvlKeyword("Width", toString(Samples()));
      o += PvlKeyword("Height", toString(Lines()));
      o += PvlKeyword("Bands", toString(Bands()));
      o += PvlKeyword("SamplingPercent", toString(SampPercent()));
      o += PvlKeyword("MinCount", toString(MinCount())); // Do we need this if EqInfo has this?

      // Create group for first file of overlap
      PvlGroup gX("File1");
      PvlKeyword stsX("StartSample", toString(StartSampleX()));
      PvlKeyword ensX("EndSample", toString(EndSampleX()));
      PvlKeyword stlX("StartLine", toString(StartLineX()));
      PvlKeyword enlX("EndLine", toString(EndLineX()));
      PvlKeyword avgX("Average");
      PvlKeyword stdX("StandardDeviation");
      PvlKeyword varX("Variance");
      for (int band = 1; band <= Bands(); band++) {
        if (HasOverlap(band)) {
          avgX += toString(GetMStats(band).X().Average());
          stdX += toString(GetMStats(band).X().StandardDeviation());
          varX += toString(GetMStats(band).X().Variance());
        }
      }
      gX += stsX;
      gX += ensX;
      gX += stlX;
      gX += enlX;
      gX += avgX;
      gX += stdX;
      gX += varX;

      // Create group for second file of overlap
      PvlGroup gY("File2");
      PvlKeyword stsY("StartSample", toString(StartSampleY()));
      PvlKeyword ensY("EndSample", toString(EndSampleY()));
      PvlKeyword stlY("StartLine", toString(StartLineY()));
      PvlKeyword enlY("EndLine", toString(EndLineY()));
      PvlKeyword avgY("Average");
      PvlKeyword stdY("StandardDeviation");
      PvlKeyword varY("Variance");
      for (int band = 1; band <= Bands(); band++) {
        if (HasOverlap(band)) {
          avgY += toString(GetMStats(band).Y().Average());
          stdY += toString(GetMStats(band).Y().StandardDeviation());
          varY += toString(GetMStats(band).Y().Variance());
        }
      }
      gY += stsY;
      gY += ensY;
      gY += stlY;
      gY += enlY;
      gY += avgY;
      gY += stdY;
      gY += varY;

      o.addGroup(gX);
      o.addGroup(gY);

      bool isValid = false;
      // Unserialize the MultivariateStatistics
      for (int band = 1; band <= Bands(); band++) {
        PvlKeyword validBand("ValidOverlap", "false");

        if (HasOverlap(band)) {
          if (IsValid(band)) {
            validBand.setValue("true");
            isValid = true;
          }
        }

        QString mStatsName = "MultivariateStatistics" + toString(band);
        PvlObject mStats(GetMStats(band).toPvl(mStatsName));
        mStats += validBand;
        o.addObject(mStats);
      }
      PvlKeyword valid("Valid", (isValid) ? "true" : "false");
      o += valid;
      return o;
    }

    catch (IException &e) {
      QString msg = "Trivial overlap between [" + FileNameX().name();
      msg += "] and [" + FileNameY().name() + "]";
      throw IException(IException::User, msg, _FILEINFO_);
    }
  }


  void OverlapStatistics::fromPvl(const PvlObject &inStats) {
    init();

    const PvlGroup &fileX = inStats.findGroup("File1");
    const PvlGroup &fileY = inStats.findGroup("File2");
    p_xFile = inStats["File1"][0];
    p_yFile = inStats["File2"][0];
    p_sampRange = inStats["Width"];
    p_lineRange = inStats["Height"];
    p_bands = inStats["Bands"];
    p_sampPercent = inStats["SamplingPercent"];

    p_minSampX = fileX["StartSample"];
    p_maxSampX = fileX["EndSample"];
    p_minLineX = fileX["StartLine"];
    p_maxLineX = fileX["EndLine"];

    p_minSampY = fileY["StartSample"];
    p_maxSampY = fileY["EndSample"];
    p_minLineY = fileY["StartLine"];
    p_maxLineY = fileY["EndLine"];

    p_mincnt = inStats["MinCount"];

    // unserialize the MStats
    for (int band = 1; band <= Bands(); band++) {
      QString name = "MultivariateStatistics" + toString(band);
      const PvlObject &mStats = inStats.findObject(name);
      p_stats.push_back(MultivariateStatistics(mStats)); 
    }

  }


  void OverlapStatistics::init() {
    p_sampPercent = 0.0;
    p_bands = 0;
    p_sampRange = 0;
    p_lineRange = 0;
    p_minSampX = 0;
    p_maxSampX = 0;
    p_minSampY = 0;
    p_maxSampY = 0;
    p_minLineX = 0;
    p_maxLineX = 0;
    p_minLineY = 0;
    p_maxLineY = 0;
    p_mincnt = 0;
  }


  std::ostream &operator<<(std::ostream &os, Isis::OverlapStatistics &stats) {
    PvlObject p = stats.toPvl();
    os << p << endl;
    return os;
  }

}