ControlNetStatistics.cpp

 
/* SPDX-License-Identifier: CC0-1.0 */
 
#include "ControlNetStatistics.h"
 
#include <QDebug>
 
#include <geos_c.h>
#include <geos/algorithm/ConvexHull.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Polygon.h>
 
#include "ControlNet.h"
#include "ControlPoint.h"
#include "ControlMeasure.h"
#include "ControlMeasureLogData.h"
#include "Cube.h"
#include "CubeManager.h"
#include "FileName.h"
#include "IString.h"
#include "Progress.h"
#include "Pvl.h"
#include "SpecialPixel.h"
#include "Statistics.h"
 
namespace Isis {
 
  QString sPointType [] = { "Fixed", "Constrained", "Free" };
 
  QString sBoolean[]    = { "False", "True" };
 
  ControlNetStatistics::ControlNetStatistics(ControlNet *pCNet, const QString &psSerialNumFile,
                                             Progress *pProgress) {
    numCNetImages = 0;
    mCNet = pCNet;
 
    mSerialNumList = SerialNumberList(psSerialNumFile);
    InitSerialNumMap();
 
    mProgress = pProgress;
 
    GetPointIntStats();
    GetPointDoubleStats();
    GenerateImageStats();
  }
 
  ControlNetStatistics::ControlNetStatistics(ControlNet *pCNet, Progress *pProgress) {
    mCNet = pCNet;
    mProgress = pProgress;
 
    GetPointIntStats();
    GetPointDoubleStats();
  }
 
  ControlNetStatistics::~ControlNetStatistics() {
    mCNet = NULL;
  }
 
  void ControlNetStatistics::InitSerialNumMap() {
    int numSn = mSerialNumList.size();
    numCNetImages = 0;
    for (int i=0; i<numSn; i++) {
      QString sn = mSerialNumList.serialNumber(i);
      mSerialNumMap[sn] = false;
    }
  }
 
  void ControlNetStatistics::GenerateControlNetStats(PvlGroup &pStatsGrp) {
    pStatsGrp = PvlGroup("ControlNetSummary");
    int numSN = mSerialNumList.size();
 
    if (numSN) {
      pStatsGrp += PvlKeyword("TotalImages",             toString(numSN));
      pStatsGrp += PvlKeyword("ImagesInControlNet", toString(numCNetImages));
    }
 
    pStatsGrp += PvlKeyword("TotalPoints",       toString(mCNet->GetNumPoints()));
    pStatsGrp += PvlKeyword("ValidPoints",       toString(NumValidPoints()));
    pStatsGrp += PvlKeyword("IgnoredPoints",     toString(mCNet->GetNumPoints() - NumValidPoints()));
    pStatsGrp += PvlKeyword("FixedPoints",       toString(NumFixedPoints()));
    pStatsGrp += PvlKeyword("ConstrainedPoints", toString(NumConstrainedPoints()));
    pStatsGrp += PvlKeyword("FreePoints",        toString(NumFreePoints()));
    pStatsGrp += PvlKeyword("EditLockPoints",    toString(mCNet->GetNumEditLockPoints()));
 
    pStatsGrp += PvlKeyword("TotalMeasures",     toString(NumMeasures()));
    pStatsGrp += PvlKeyword("ValidMeasures",     toString(NumValidMeasures()));
    pStatsGrp += PvlKeyword("IgnoredMeasures",   toString(NumIgnoredMeasures()));
    pStatsGrp += PvlKeyword("EditLockMeasures",  toString(mCNet->GetNumEditLockMeasures()));
 
    double dValue = GetAverageResidual();
    pStatsGrp += PvlKeyword("AvgResidual",       (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMinimumResidual();
    pStatsGrp += PvlKeyword("MinResidual",       (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMaximumResidual();
    pStatsGrp += PvlKeyword("MaxResidual",       (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMinLineResidual();
    pStatsGrp += PvlKeyword("MinLineResidual",   (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMaxLineResidual();
    pStatsGrp += PvlKeyword("MaxLineResidual",   (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMinSampleResidual();
    pStatsGrp += PvlKeyword("MinSampleResidual", (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMaxSampleResidual();
    pStatsGrp += PvlKeyword("MaxSampleResidual", (dValue == Null ? "Null" : toString(dValue)));
 
    // Shifts - Line, Sample, Pixel
    dValue = GetMinLineShift();
    pStatsGrp += PvlKeyword("MinLineShift",      (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMaxLineShift();
    pStatsGrp += PvlKeyword("MaxLineShift",      (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMinSampleShift();
    pStatsGrp += PvlKeyword("MinSampleShift",    (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetMaxSampleShift();
    pStatsGrp += PvlKeyword("MaxSampleShift",    (dValue == Null ? "Null" : toString(dValue)));
 
    dValue = GetAvgPixelShift();
    pStatsGrp += PvlKeyword("AvgPixelShift",     (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = GetMinPixelShift();
    pStatsGrp += PvlKeyword("MinPixelShift",     (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = GetMaxPixelShift();
    pStatsGrp += PvlKeyword("MaxPixelShift",     (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = mPointDoubleStats[minGFit];
    pStatsGrp += PvlKeyword("MinGoodnessOfFit",  (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = mPointDoubleStats[maxGFit];
    pStatsGrp += PvlKeyword("MaxGoodnessOfFit",  (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = mPointDoubleStats[minEccentricity];
    pStatsGrp += PvlKeyword("MinEccentricity",   (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = mPointDoubleStats[maxEccentricity];
    pStatsGrp += PvlKeyword("MaxEccentricity",   (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = mPointDoubleStats[minPixelZScore];
    pStatsGrp += PvlKeyword("MinPixelZScore",    (dValue == Null ? "NA" : toString(dValue)));
 
    dValue = mPointDoubleStats[maxPixelZScore];
    pStatsGrp += PvlKeyword("MaxPixelZScore",    (dValue == Null ? "NA" : toString(dValue)));
 
    // Convex Hull
    if (mSerialNumList.size()) {
      dValue = mConvexHullRatioStats.Minimum();
      pStatsGrp += PvlKeyword("MinConvexHullRatio", (dValue == Null ? "Null" : toString(dValue)));
 
      dValue = mConvexHullRatioStats.Maximum();
      pStatsGrp += PvlKeyword("MaxConvexHullRatio", (dValue == Null ? "Null" : toString(dValue)));
 
      dValue = mConvexHullRatioStats.Average();
      pStatsGrp += PvlKeyword("AvgConvexHullRatio", (dValue == Null ? "Null" : toString(dValue)));
    }
  }
 
 
  void ControlNetStatistics::GenerateImageStats() {
    geos::geom::GeometryFactory::Ptr geosFactory = geos::geom::GeometryFactory::create();
 
    CubeManager cubeMgr;
    cubeMgr.SetNumOpenCubes(50);
 
    QList<QString> cnetSerials = mCNet->GetCubeSerials();
 
    if (mProgress != NULL) {
      mProgress->SetText("Generating Image Stats.....");
      mProgress->SetMaximumSteps(cnetSerials.size());
      mProgress->CheckStatus();
    }
 
    foreach (QString sn, cnetSerials) {
      geos::geom::CoordinateSequence ptCoordinates;
 
      // setup vector for number of image properties and init to 0
      QVector<double> imgStats(numImageStats, 0);
 
      // Open the cube to get the dimensions
      Cube *cube = cubeMgr.OpenCube(mSerialNumList.fileName(sn));
 
      mSerialNumMap[sn] = true;
      numCNetImages++;
 
      imgStats[imgSamples] = cube->sampleCount();
      imgStats[imgLines]   = cube->lineCount();
      double cubeArea      = imgStats[imgSamples] * imgStats[imgLines];
 
      QList< ControlMeasure * > measures = mCNet->GetMeasuresInCube(sn);
 
      // Populate pts with a list of control points
      if (!measures.isEmpty()) {
        foreach (ControlMeasure * measure, measures) {
          ControlPoint *parentPoint = measure->Parent();
          imgStats[imgTotalPoints]++;
          if (parentPoint->IsIgnored()) {
            imgStats[imgIgnoredPoints]++;
          }
          if (parentPoint->GetType() == ControlPoint::Fixed) {
            imgStats[imgFixedPoints]++;
          }
          if (parentPoint->GetType() == ControlPoint::Constrained) {
            imgStats[imgConstrainedPoints]++;
          }
          if (parentPoint->GetType() == ControlPoint::Free) {
            imgStats[imgFreePoints]++;
          }
          if (parentPoint->IsEditLocked()) {
            imgStats[imgLockedPoints]++;
          }
          if (measure->IsEditLocked()) {
            imgStats[imgLocked]++;
          }
          ptCoordinates.add(geos::geom::Coordinate(measure->GetSample(),
                                                    measure->GetLine()));
        }
 
        ptCoordinates.add(geos::geom::Coordinate(measures[0]->GetSample(),
                                                  measures[0]->GetLine()));
      }
 
      if (ptCoordinates.size() >= 4) {
        // Calculate the convex hull
 
        // Even though geos doesn't create valid linear rings/polygons from this set of coordinates,
        //   because it self-intersects many many times, it still correctly does a convex hull
        //   calculation on the points in the polygon.
        geos::geom::Geometry * convexHull = geosFactory->createPolygon(
          geosFactory->createLinearRing(ptCoordinates))->convexHull().release();
 
        // Calculate the area of the convex hull
        imgStats[imgConvexHullArea] = convexHull->getArea();
        imgStats[imgConvexHullRatio] = imgStats[imgConvexHullArea] / cubeArea;
      }
 
      // Add info to statistics to get min, max and avg convex hull
      mConvexHullStats.AddData(imgStats[imgConvexHullArea]);
      mConvexHullRatioStats.AddData(imgStats[imgConvexHullRatio]);
 
      mImageMap[sn] = imgStats;
 
      // Update Progress
      if (mProgress != NULL)
        mProgress->CheckStatus();
    }
  }
 
 
  void ControlNetStatistics::PrintImageStats(const QString &psImageFile) {
    // Check if the image list has been provided
    if (!mSerialNumList.size()) {
      QString msg = "Serial Number of Images has not been provided to get Image Stats";
      throw IException(IException::User, msg, _FILEINFO_);
    }
 
    FileName outFile(psImageFile);
    ofstream ostm;
    QString outName(outFile.expanded());
    ostm.open(outName.toLatin1().data(), std::ios::out);
 
    if ( ostm.fail() ) {
      QString msg = QObject::tr("Cannot open file [%1]").arg(psImageFile);
      throw IException(IException::Io, msg, _FILEINFO_);
     }
 
    //map< QString, vector<double> >::iterator it;
    // imgSamples, imgLines, imgTotalPoints, imgIgnoredPoints, imgFixedPoints, imgLockedPoints, imgLocked, imgConstrainedPoints, imgFreePoints, imgConvexHullArea, imgConvexHullRatio
    QMap<QString, bool>::iterator it;
    // Log into the output file
    ostm << "Filename, SerialNumber, TotalPoints, PointsIgnored, PointsEditLocked, Fixed, Constrained, Free, ConvexHullRatio" <<  endl;
    //for (it = mImageMap.begin(); it != mImageMap.end(); it++) {
 
    for (it = mSerialNumMap.begin(); it != mSerialNumMap.end(); it++) {
      ostm << mSerialNumList.fileName(it.key()) << ", " << it.key() << ", ";
      bool serialNumExists = it.value();
      if (serialNumExists) {
        QVector<double> imgStats = mImageMap[(it).key()] ;
        ostm << imgStats[imgTotalPoints]<< ", " << imgStats[imgIgnoredPoints] << ", " ;
        ostm << imgStats[imgLockedPoints] << ", " << imgStats[imgFixedPoints] << ", " ;
        ostm << imgStats[imgConstrainedPoints] << ", " << imgStats[imgFreePoints] << ", ";
        ostm << imgStats[imgConvexHullRatio] << endl;
      }
      else {
        ostm << "0, 0, 0, 0, 0, 0, 0" << endl;
      }
    }
 
    if (!ostm) {
      QString msg = QObject::tr("Error writing to file: [%1]").arg(psImageFile);
      throw IException(IException::Io, msg, _FILEINFO_);
    }
    ostm.close();
  }
 
 
  QVector<double> ControlNetStatistics::GetImageStatsBySerialNum(QString psSerialNum) const {
    return (mImageMap.find(psSerialNum)).value();
  }
 
 
  void ControlNetStatistics::GeneratePointStats(const QString &psPointFile) {
    Isis::FileName outFile(psPointFile);
 
    ofstream ostm;
    QString outName(outFile.expanded());
    ostm.open(outName.toLatin1().data(), std::ios::out);
 
    if ( ostm.fail() ) {
      QString msg = QObject::tr("Cannot open file [%1]").arg(psPointFile);
      throw IException(IException::Io, msg, _FILEINFO_);
    }
 
    ostm << " PointId, PointType, PointIgnore, PointEditLock, TotalMeasures, MeasuresValid, MeasuresIgnore, MeasuresEditLock" << endl;
 
    int iNumPoints = mCNet->GetNumPoints();
 
    // Initialise the Progress object
    if (mProgress != NULL && iNumPoints > 0) {
      mProgress->SetText("Point Stats: Loading Control Points...");
      mProgress->SetMaximumSteps(iNumPoints);
      mProgress->CheckStatus();
    }
 
    for (int i = 0; i < iNumPoints; i++) {
      const ControlPoint *cPoint = mCNet->GetPoint(i);
      int iNumMeasures     = cPoint->GetNumMeasures();
      int iValidMeasures   = cPoint->GetNumValidMeasures();
      int iIgnoredMeasures = iNumMeasures - iValidMeasures;
 
      // Log into the output file
      ostm << cPoint->GetId()   << ", " << sPointType[(int)cPoint->GetType()] << ", " << sBoolean[(int)cPoint->IsIgnored()] << ", " ;
      ostm << sBoolean[(int)cPoint->IsEditLocked()] << ", " << iNumMeasures << ", " << iValidMeasures << ", ";
      ostm << iIgnoredMeasures << ", " << cPoint->GetNumLockedMeasures() << endl;
 
      // Update Progress
      if (mProgress != NULL)
        mProgress->CheckStatus();
    }
 
    if (!ostm) {
      QString msg = QObject::tr("Error writing to file: [%1]").arg(psPointFile);
      throw IException(IException::Io, msg, _FILEINFO_);
    }
 
    ostm.close();
  }
 
 
  void ControlNetStatistics::GetPointIntStats() {
    // Init all the entries
    // totalPoints, validPoints, ignoredPoints, fixedPoints, constrainedPoints, editLockedPoints,
    // totalMeasures, validMeasures, ignoredMeasures, editLockedMeasures
    for (int i=0; i<numPointIntStats; i++) {
      mPointIntStats[i] = 0;
    }
 
    int iNumPoints = mCNet->GetNumPoints();
 
    // totalPoints
    mPointIntStats[totalPoints] = iNumPoints;
 
    for (int i = 0; i < iNumPoints; i++) {
      if (!mCNet->GetPoint(i)->IsIgnored()) {
        // validPoints
        mPointIntStats[validPoints]++;
      }
      else {
        // ignoredPoints
        mPointIntStats[ignoredPoints]++;
      }
 
      // fixedPoints
      if (mCNet->GetPoint(i)->GetType() == ControlPoint::Fixed)
        mPointIntStats[fixedPoints]++;
 
      // constrainedPoints
      if (mCNet->GetPoint(i)->GetType() == ControlPoint::Constrained)
        mPointIntStats[constrainedPoints]++;
 
      // free points
      if (mCNet->GetPoint(i)->GetType() == ControlPoint::Free)
        mPointIntStats[freePoints]++;
 
      // editLockedPoints
      if (mCNet->GetPoint(i)->IsEditLocked()) {
        mPointIntStats[editLockedPoints]++;
      }
 
      // totalMeasures
      mPointIntStats[totalMeasures] += mCNet->GetPoint(i)->GetNumMeasures();
 
      // validMeasures
      mPointIntStats[validMeasures] += mCNet->GetPoint(i)->GetNumValidMeasures();
 
      // editLockedMeasures
      mPointIntStats[editLockedMeasures] +=  mCNet->GetPoint(i)->GetNumLockedMeasures();
    }
    // ignoredMeasures
    mPointIntStats[ignoredMeasures] = mPointIntStats[totalMeasures] -  mPointIntStats[validMeasures];
  }
 
 
  void ControlNetStatistics::InitPointDoubleStats() {
    for (int i = 0; i < numPointDblStats; i++) {
      mPointDoubleStats[i] = Null;
    }
  }
 
 
  void ControlNetStatistics::GetPointDoubleStats() {
    InitPointDoubleStats();
 
    int iNumPoints = mCNet->GetNumPoints();
    double dValue = 0;
 
    Statistics residualMagStats;
    Statistics pixelShiftStats;
 
    for (int i = 0; i < iNumPoints; i++) {
      ControlPoint * cp = mCNet->GetPoint(i);
 
      if (!cp->IsIgnored()) {
        for (int cmIndex = 0; cmIndex < cp->GetNumMeasures(); cmIndex++) {
          ControlMeasure *cm = cp->GetMeasure(cmIndex);
 
          if (!cm->IsIgnored()) {
            residualMagStats.AddData(cm->GetResidualMagnitude());
 
            if (!IsSpecial(cm->GetPixelShift()))
              pixelShiftStats.AddData(fabs(cm->GetPixelShift()));
          }
        }
      }
 
      Statistics resMagStats = cp->GetStatistic(
          &ControlMeasure::GetResidualMagnitude);
      UpdateMinMaxStats(resMagStats, minResidual, maxResidual);
 
      Statistics resLineStats = cp->GetStatistic(
          &ControlMeasure::GetLineResidual);
      UpdateMinMaxStats(resLineStats, minLineResidual, maxLineResidual);
 
      Statistics resSampStats = cp->GetStatistic(
          &ControlMeasure::GetSampleResidual);
      UpdateMinMaxStats(resSampStats, minSampleResidual, maxSampleResidual);
 
      Statistics pixShiftStats = cp->GetStatistic(
          &ControlMeasure::GetPixelShift);
      UpdateMinMaxStats(pixShiftStats, minPixelShift, maxPixelShift);
 
      Statistics lineShiftStats = cp->GetStatistic(
          &ControlMeasure::GetLineShift);
      UpdateMinMaxStats(lineShiftStats, minLineShift, maxLineShift);
 
      Statistics sampShiftStats = cp->GetStatistic(
          &ControlMeasure::GetSampleShift);
      UpdateMinMaxStats(sampShiftStats, minSampleShift, maxSampleShift);
 
      Statistics gFitStats = cp->GetStatistic(
          ControlMeasureLogData::GoodnessOfFit);
      UpdateMinMaxStats(gFitStats, minGFit, maxGFit);
 
      Statistics minPixelZScoreStats = cp->GetStatistic(
          ControlMeasureLogData::MinimumPixelZScore);
 
      if (minPixelZScoreStats.ValidPixels()) {
        dValue = fabs(minPixelZScoreStats.Minimum());
        if (mPointDoubleStats[minPixelZScore] > dValue)
          mPointDoubleStats[minPixelZScore] = dValue;
      }
 
      Statistics maxPixelZScoreStats = cp->GetStatistic(
          ControlMeasureLogData::MaximumPixelZScore);
 
      if (maxPixelZScoreStats.ValidPixels()) {
        dValue = fabs(maxPixelZScoreStats.Maximum());
        if (mPointDoubleStats[maxPixelZScore] > dValue)
          mPointDoubleStats[maxPixelZScore] = dValue;
      }
    }
 
    // Average Residuals
    mPointDoubleStats[avgResidual] = residualMagStats.Average();
 
    // Average Shift
    mPointDoubleStats[avgPixelShift] = pixelShiftStats.Average();
  }
 
 
  void ControlNetStatistics::UpdateMinMaxStats(const Statistics & stats,
      ePointDoubleStats min, ePointDoubleStats max) {
    if (stats.ValidPixels()) {
      if (mPointDoubleStats[min] != Null) {
        mPointDoubleStats[min] = qMin(
            mPointDoubleStats[min], fabs(stats.Minimum()));
      }
      else {
        mPointDoubleStats[min] = fabs(stats.Minimum());
      }
 
      if (mPointDoubleStats[max] != Null) {
        mPointDoubleStats[max] = qMax(
            mPointDoubleStats[max], fabs(stats.Maximum()));
      }
      else {
        mPointDoubleStats[max] = fabs(stats.Maximum());
      }
    }
  }
}