USGS

Isis 3.0 Object Programmers' Reference

 Home
ControlNetStatistics.cpp
1 #include "ControlNetStatistics.h"
2 
3 #include <QDebug>
4 
5 #include <geos_c.h>
6 #include <geos/algorithm/ConvexHull.h>
7 #include <geos/geom/CoordinateSequence.h>
8 #include <geos/geom/CoordinateArraySequence.h>
9 #include <geos/geom/Envelope.h>
10 #include <geos/geom/Geometry.h>
11 #include <geos/geom/GeometryFactory.h>
12 #include <geos/geom/Polygon.h>
13 
14 #include "ControlNet.h"
15 #include "ControlPoint.h"
16 #include "ControlMeasure.h"
17 #include "ControlMeasureLogData.h"
18 #include "ControlCubeGraphNode.h"
19 #include "Cube.h"
20 #include "CubeManager.h"
21 #include "FileName.h"
22 #include "IString.h"
23 #include "Progress.h"
24 #include "Pvl.h"
25 #include "SpecialPixel.h"
26 #include "Statistics.h"
27 
28 
29 using namespace std;
30 
31 namespace Isis {
32 
34  QString sPointType [] = { "Fixed", "Constrained", "Free" };
35 
37  QString sBoolean[] = { "False", "True" };
38 
48  ControlNetStatistics::ControlNetStatistics(ControlNet *pCNet, const QString &psSerialNumFile,
49  Progress *pProgress) {
50  numCNetImages = 0;
51  mCNet = pCNet;
52 
53  mSerialNumList = SerialNumberList(psSerialNumFile);
54  InitSerialNumMap();
55 
56  mProgress = pProgress;
57 
58  GetPointIntStats();
59  GetPointDoubleStats();
60  GenerateImageStats();
61  }
62 
71  ControlNetStatistics::ControlNetStatistics(ControlNet *pCNet, Progress *pProgress) {
72  mCNet = pCNet;
73  mProgress = pProgress;
74 
75  GetPointIntStats();
76  GetPointDoubleStats();
77  }
78 
84  ControlNetStatistics::~ControlNetStatistics() {
85  mCNet = NULL;
86  }
87 
93  void ControlNetStatistics::InitSerialNumMap() {
94  int numSn = mSerialNumList.size();
95  numCNetImages = 0;
96  for (int i=0; i<numSn; i++) {
97  QString sn = mSerialNumList.serialNumber(i);
98  mSerialNumMap[sn] = false;
99  }
100  }
101 
112  void ControlNetStatistics::GenerateControlNetStats(PvlGroup &pStatsGrp) {
113  pStatsGrp = PvlGroup("ControlNetSummary");
114  int numSN = mSerialNumList.size();
115 
116  if (numSN) {
117  pStatsGrp += PvlKeyword("TotalImages", toString(numSN));
118  pStatsGrp += PvlKeyword("ImagesInControlNet", toString(numCNetImages));
119  }
120 
121  pStatsGrp += PvlKeyword("TotalPoints", toString(mCNet->GetNumPoints()));
122  pStatsGrp += PvlKeyword("ValidPoints", toString(NumValidPoints()));
123  pStatsGrp += PvlKeyword("IgnoredPoints", toString(mCNet->GetNumPoints() - NumValidPoints()));
124  pStatsGrp += PvlKeyword("FixedPoints", toString(NumFixedPoints()));
125  pStatsGrp += PvlKeyword("ConstrainedPoints", toString(NumConstrainedPoints()));
126  pStatsGrp += PvlKeyword("FreePoints", toString(NumFreePoints()));
127  pStatsGrp += PvlKeyword("EditLockPoints", toString(mCNet->GetNumEditLockPoints()));
128 
129  pStatsGrp += PvlKeyword("TotalMeasures", toString(NumMeasures()));
130  pStatsGrp += PvlKeyword("ValidMeasures", toString(NumValidMeasures()));
131  pStatsGrp += PvlKeyword("IgnoredMeasures", toString(NumIgnoredMeasures()));
132  pStatsGrp += PvlKeyword("EditLockMeasures", toString(mCNet->GetNumEditLockMeasures()));
133 
134  double dValue = GetAverageResidual();
135  pStatsGrp += PvlKeyword("AvgResidual", (dValue == Null ? "Null" : toString(dValue)));
136 
137  dValue = GetMinimumResidual();
138  pStatsGrp += PvlKeyword("MinResidual", (dValue == Null ? "Null" : toString(dValue)));
139 
140  dValue = GetMaximumResidual();
141  pStatsGrp += PvlKeyword("MaxResidual", (dValue == Null ? "Null" : toString(dValue)));
142 
143  dValue = GetMinLineResidual();
144  pStatsGrp += PvlKeyword("MinLineResidual", (dValue == Null ? "Null" : toString(dValue)));
145 
146  dValue = GetMaxLineResidual();
147  pStatsGrp += PvlKeyword("MaxLineResidual", (dValue == Null ? "Null" : toString(dValue)));
148 
149  dValue = GetMinSampleResidual();
150  pStatsGrp += PvlKeyword("MinSampleResidual", (dValue == Null ? "Null" : toString(dValue)));
151 
152  dValue = GetMaxSampleResidual();
153  pStatsGrp += PvlKeyword("MaxSampleResidual", (dValue == Null ? "Null" : toString(dValue)));
154 
155  // Shifts - Line, Sample, Pixel
156  dValue = GetMinLineShift();
157  pStatsGrp += PvlKeyword("MinLineShift", (dValue == Null ? "Null" : toString(dValue)));
158 
159  dValue = GetMaxLineShift();
160  pStatsGrp += PvlKeyword("MaxLineShift", (dValue == Null ? "Null" : toString(dValue)));
161 
162  dValue = GetMinSampleShift();
163  pStatsGrp += PvlKeyword("MinSampleShift", (dValue == Null ? "Null" : toString(dValue)));
164 
165  dValue = GetMaxSampleShift();
166  pStatsGrp += PvlKeyword("MaxSampleShift", (dValue == Null ? "Null" : toString(dValue)));
167 
168  dValue = GetAvgPixelShift();
169  pStatsGrp += PvlKeyword("AvgPixelShift", (dValue == Null ? "NA" : toString(dValue)));
170 
171  dValue = GetMinPixelShift();
172  pStatsGrp += PvlKeyword("MinPixelShift", (dValue == Null ? "NA" : toString(dValue)));
173 
174  dValue = GetMaxPixelShift();
175  pStatsGrp += PvlKeyword("MaxPixelShift", (dValue == Null ? "NA" : toString(dValue)));
176 
177  dValue = mPointDoubleStats[minGFit];
178  pStatsGrp += PvlKeyword("MinGoodnessOfFit", (dValue == Null ? "NA" : toString(dValue)));
179 
180  dValue = mPointDoubleStats[maxGFit];
181  pStatsGrp += PvlKeyword("MaxGoodnessOfFit", (dValue == Null ? "NA" : toString(dValue)));
182 
183  dValue = mPointDoubleStats[minEccentricity];
184  pStatsGrp += PvlKeyword("MinEccentricity", (dValue == Null ? "NA" : toString(dValue)));
185 
186  dValue = mPointDoubleStats[maxEccentricity];
187  pStatsGrp += PvlKeyword("MaxEccentricity", (dValue == Null ? "NA" : toString(dValue)));
188 
189  dValue = mPointDoubleStats[minPixelZScore];
190  pStatsGrp += PvlKeyword("MinPixelZScore", (dValue == Null ? "NA" : toString(dValue)));
191 
192  dValue = mPointDoubleStats[maxPixelZScore];
193  pStatsGrp += PvlKeyword("MaxPixelZScore", (dValue == Null ? "NA" : toString(dValue)));
194 
195  // Convex Hull
196  if (mSerialNumList.size()) {
197  dValue = mConvexHullRatioStats.Minimum();
198  pStatsGrp += PvlKeyword("MinConvexHullRatio", (dValue == Null ? "Null" : toString(dValue)));
199 
200  dValue = mConvexHullRatioStats.Maximum();
201  pStatsGrp += PvlKeyword("MaxConvexHullRatio", (dValue == Null ? "Null" : toString(dValue)));
202 
203  dValue = mConvexHullRatioStats.Average();
204  pStatsGrp += PvlKeyword("AvgConvexHullRatio", (dValue == Null ? "Null" : toString(dValue)));
205  }
206  }
207 
208 
216  void ControlNetStatistics::GenerateImageStats() {
217  geos::geom::GeometryFactory geosFactory;
218 
219  CubeManager cubeMgr;
220  cubeMgr.SetNumOpenCubes(50);
221 
222  mCubeGraphNodes = mCNet->GetCubeGraphNodes();
223 
224  if (mProgress != NULL) {
225  mProgress->SetText("Generating Image Stats.....");
226  mProgress->SetMaximumSteps(mCubeGraphNodes.size());
227  mProgress->CheckStatus();
228  }
229 
230  foreach (ControlCubeGraphNode * node, mCubeGraphNodes) {
231  geos::geom::CoordinateSequence * ptCoordinates =
232  new geos::geom::CoordinateArraySequence();
233 
234  // setup vector for number of image properties and init to 0
235  vector<double> imgStats(numImageStats, 0);
236 
237  // Open the cube to get the dimensions
238  QString sn = node->getSerialNumber();
239  Cube *cube = cubeMgr.OpenCube(mSerialNumList.fileName(sn));
240 
241  mSerialNumMap[sn] = true;
242  numCNetImages++;
243 
244  imgStats[imgSamples] = cube->sampleCount();
245  imgStats[imgLines] = cube->lineCount();
246  double cubeArea = imgStats[imgSamples] * imgStats[imgLines];
247 
248  QList< ControlMeasure * > measures = node->getMeasures();
249 
250  // Populate pts with a list of control points
251  if (!measures.isEmpty()) {
252  foreach (ControlMeasure * measure, measures) {
253  ControlPoint *parentPoint = measure->Parent();
254  imgStats[imgTotalPoints]++;
255  if (parentPoint->IsIgnored()) {
256  imgStats[imgIgnoredPoints]++;
257  }
258  if (parentPoint->GetType() == ControlPoint::Fixed) {
259  imgStats[imgFixedPoints]++;
260  }
261  if (parentPoint->GetType() == ControlPoint::Constrained) {
262  imgStats[imgConstrainedPoints]++;
263  }
264  if (parentPoint->GetType() == ControlPoint::Free) {
265  imgStats[imgFreePoints]++;
266  }
267  if (parentPoint->IsEditLocked()) {
268  imgStats[imgLockedPoints]++;
269  }
270  if (measure->IsEditLocked()) {
271  imgStats[imgLocked]++;
272  }
273  ptCoordinates->add(geos::geom::Coordinate(measure->GetSample(),
274  measure->GetLine()));
275  }
276 
277  ptCoordinates->add(geos::geom::Coordinate(measures[0]->GetSample(),
278  measures[0]->GetLine()));
279  }
280 
281  if (ptCoordinates->size() >= 4) {
282  // Calculate the convex hull
283 
284  // Even though geos doesn't create valid linear rings/polygons from this set of coordinates,
285  // because it self-intersects many many times, it still correctly does a convex hull
286  // calculation on the points in the polygon.
287  geos::geom::Geometry * convexHull = geosFactory.createPolygon(
288  geosFactory.createLinearRing(ptCoordinates), 0)->convexHull();
289 
290  // Calculate the area of the convex hull
291  imgStats[imgConvexHullArea] = convexHull->getArea();
292  imgStats[imgConvexHullRatio] = imgStats[imgConvexHullArea] / cubeArea;
293  }
294 
295  // Add info to statistics to get min, max and avg convex hull
296  mConvexHullStats.AddData(imgStats[imgConvexHullArea]);
297  mConvexHullRatioStats.AddData(imgStats[imgConvexHullRatio]);
298 
299  mImageMap[sn] = imgStats;
300 
301  delete ptCoordinates;
302  ptCoordinates = NULL;
303 
304  // Update Progress
305  if (mProgress != NULL)
306  mProgress->CheckStatus();
307  }
308  }
309 
310 
320  void ControlNetStatistics::PrintImageStats(const QString &psImageFile) {
321  // Check if the image list has been provided
322  if (!mSerialNumList.size()) {
323  QString msg = "Serial Number of Images has not been provided to get Image Stats";
324  throw IException(IException::User, msg, _FILEINFO_);
325  }
326 
327  FileName outFile(psImageFile);
328  ofstream ostm;
329  QString outName(outFile.expanded());
330  ostm.open(outName.toAscii().data(), std::ios::out);
331 
332  if ( ostm.fail() ) {
333  QString msg = QObject::tr("Cannot open file [%1]").arg(psImageFile);
334  throw IException(IException::Io, msg, _FILEINFO_);
335  }
336 
337  //map< QString, vector<double> >::iterator it;
338  map<QString, bool>::iterator it;
339  // imgSamples, imgLines, imgTotalPoints, imgIgnoredPoints, imgFixedPoints, imgLockedPoints, imgLocked, imgConstrainedPoints, imgFreePoints, imgConvexHullArea, imgConvexHullRatio
340 
341  // Log into the output file
342  ostm << "Filename, SerialNumber, TotalPoints, PointsIgnored, PointsEditLocked, Fixed, Constrained, Free, ConvexHullRatio" << endl;
343  //for (it = mImageMap.begin(); it != mImageMap.end(); it++) {
344  for (it = mSerialNumMap.begin(); it != mSerialNumMap.end(); it++) {
345  ostm << mSerialNumList.fileName((*it).first) << ", " << (*it).first << ", ";
346  bool serialNumExists = (*it).second ;
347  if (serialNumExists) {
348  vector<double>imgStats = mImageMap[(*it).first] ;
349  ostm << imgStats[imgTotalPoints]<< ", " << imgStats[imgIgnoredPoints] << ", " ;
350  ostm << imgStats[imgLockedPoints] << ", " << imgStats[imgFixedPoints] << ", " ;
351  ostm << imgStats[imgConstrainedPoints] << ", " << imgStats[imgFreePoints] << ", ";
352  ostm << imgStats[ imgConvexHullRatio] << endl;
353  }
354  else {
355  ostm << "0, 0, 0, 0, 0, 0, 0" << endl;
356  }
357  }
358 
359  if (!ostm) {
360  QString msg = QObject::tr("Error writing to file: [%1]").arg(psImageFile);
361  throw IException(IException::Io, msg, _FILEINFO_);
362  }
363  ostm.close();
364  }
365 
366 
376  vector<double> ControlNetStatistics::GetImageStatsBySerialNum(QString psSerialNum) const {
377  return (*mImageMap.find(psSerialNum)).second;
378  }
379 
380 
390  void ControlNetStatistics::GeneratePointStats(const QString &psPointFile) {
391  Isis::FileName outFile(psPointFile);
392 
393  ofstream ostm;
394  QString outName(outFile.expanded());
395  ostm.open(outName.toAscii().data(), std::ios::out);
396 
397  if ( ostm.fail() ) {
398  QString msg = QObject::tr("Cannot open file [%1]").arg(psPointFile);
399  throw IException(IException::Io, msg, _FILEINFO_);
400  }
401 
402  ostm << " PointId, PointType, PointIgnore, PointEditLock, TotalMeasures, MeasuresValid, MeasuresIgnore, MeasuresEditLock," << endl;
403 
404  int iNumPoints = mCNet->GetNumPoints();
405 
406  // Initialise the Progress object
407  if (mProgress != NULL && iNumPoints > 0) {
408  mProgress->SetText("Point Stats: Loading Control Points...");
409  mProgress->SetMaximumSteps(iNumPoints);
410  mProgress->CheckStatus();
411  }
412 
413  for (int i = 0; i < iNumPoints; i++) {
414  const ControlPoint *cPoint = mCNet->GetPoint(i);
415  int iNumMeasures = cPoint->GetNumMeasures();
416  int iValidMeasures = cPoint->GetNumValidMeasures();
417  int iIgnoredMeasures = iNumMeasures - iValidMeasures;
418 
419  // Log into the output file
420  ostm << cPoint->GetId() << ", " << sPointType[(int)cPoint->GetType()] << ", " << sBoolean[(int)cPoint->IsIgnored()] << ", " ;
421  ostm << sBoolean[(int)cPoint->IsEditLocked()] << ", " << iNumMeasures << ", " << iValidMeasures << ", ";
422  ostm << iIgnoredMeasures << ", " << cPoint->GetNumLockedMeasures() << endl;
423 
424  // Update Progress
425  if (mProgress != NULL)
426  mProgress->CheckStatus();
427  }
428 
429  if (!ostm) {
430  QString msg = QObject::tr("Error writing to file: [%1]").arg(psPointFile);
431  throw IException(IException::Io, msg, _FILEINFO_);
432  }
433 
434  ostm.close();
435  }
436 
437 
443  void ControlNetStatistics::GetPointIntStats() {
444  // Init all the entries
445  // totalPoints, validPoints, ignoredPoints, fixedPoints, constrainedPoints, editLockedPoints,
446  // totalMeasures, validMeasures, ignoredMeasures, editLockedMeasures
447  for (int i=0; i<numPointIntStats; i++) {
448  mPointIntStats[i] = 0;
449  }
450 
451  int iNumPoints = mCNet->GetNumPoints();
452 
453  // totalPoints
454  mPointIntStats[totalPoints] = iNumPoints;
455 
456  for (int i = 0; i < iNumPoints; i++) {
457  if (!mCNet->GetPoint(i)->IsIgnored()) {
458  // validPoints
459  mPointIntStats[validPoints]++;
460  }
461  else {
462  // ignoredPoints
463  mPointIntStats[ignoredPoints]++;
464  }
465 
466  // fixedPoints
467  if (mCNet->GetPoint(i)->GetType() == ControlPoint::Fixed)
468  mPointIntStats[fixedPoints]++;
469 
470  // constrainedPoints
471  if (mCNet->GetPoint(i)->GetType() == ControlPoint::Constrained)
472  mPointIntStats[constrainedPoints]++;
473 
474  // free points
475  if (mCNet->GetPoint(i)->GetType() == ControlPoint::Free)
476  mPointIntStats[freePoints]++;
477 
478  // editLockedPoints
479  if (mCNet->GetPoint(i)->IsEditLocked()) {
480  mPointIntStats[editLockedPoints]++;
481  }
482 
483  // totalMeasures
484  mPointIntStats[totalMeasures] += mCNet->GetPoint(i)->GetNumMeasures();
485 
486  // validMeasures
487  mPointIntStats[validMeasures] += mCNet->GetPoint(i)->GetNumValidMeasures();
488 
489  // editLockedMeasures
490  mPointIntStats[editLockedMeasures] += mCNet->GetPoint(i)->GetNumLockedMeasures();
491  }
492  // ignoredMeasures
493  mPointIntStats[ignoredMeasures] = mPointIntStats[totalMeasures] - mPointIntStats[validMeasures];
494  }
495 
496 
502  void ControlNetStatistics::InitPointDoubleStats() {
503  for (int i = 0; i < numPointDblStats; i++) {
504  mPointDoubleStats[i] = Null;
505  }
506  }
507 
508 
515  void ControlNetStatistics::GetPointDoubleStats() {
516  InitPointDoubleStats();
517 
518  int iNumPoints = mCNet->GetNumPoints();
519  double dValue = 0;
520 
521  Statistics residualMagStats;
522  Statistics pixelShiftStats;
523 
524  for (int i = 0; i < iNumPoints; i++) {
525  ControlPoint * cp = mCNet->GetPoint(i);
526 
527  if (!cp->IsIgnored()) {
528  for (int cmIndex = 0; cmIndex < cp->GetNumMeasures(); cmIndex++) {
529  ControlMeasure *cm = cp->GetMeasure(cmIndex);
530 
531  if (!cm->IsIgnored()) {
532  residualMagStats.AddData(cm->GetResidualMagnitude());
533 
534  if (!IsSpecial(cm->GetPixelShift()))
535  pixelShiftStats.AddData(fabs(cm->GetPixelShift()));
536  }
537  }
538  }
539 
540  Statistics resMagStats = cp->GetStatistic(
541  &ControlMeasure::GetResidualMagnitude);
542  UpdateMinMaxStats(resMagStats, minResidual, maxResidual);
543 
544  Statistics resLineStats = cp->GetStatistic(
545  &ControlMeasure::GetLineResidual);
546  UpdateMinMaxStats(resLineStats, minLineResidual, maxLineResidual);
547 
548  Statistics resSampStats = cp->GetStatistic(
549  &ControlMeasure::GetSampleResidual);
550  UpdateMinMaxStats(resSampStats, minSampleResidual, maxSampleResidual);
551 
552  Statistics pixShiftStats = cp->GetStatistic(
553  &ControlMeasure::GetPixelShift);
554  UpdateMinMaxStats(pixShiftStats, minPixelShift, maxPixelShift);
555 
556  Statistics lineShiftStats = cp->GetStatistic(
557  &ControlMeasure::GetLineShift);
558  UpdateMinMaxStats(lineShiftStats, minLineShift, maxLineShift);
559 
560  Statistics sampShiftStats = cp->GetStatistic(
561  &ControlMeasure::GetSampleShift);
562  UpdateMinMaxStats(sampShiftStats, minSampleShift, maxSampleShift);
563 
564  Statistics gFitStats = cp->GetStatistic(
565  ControlMeasureLogData::GoodnessOfFit);
566  UpdateMinMaxStats(gFitStats, minGFit, maxGFit);
567 
568  Statistics minPixelZScoreStats = cp->GetStatistic(
569  ControlMeasureLogData::MinimumPixelZScore);
570 
571  if (minPixelZScoreStats.ValidPixels()) {
572  dValue = fabs(minPixelZScoreStats.Minimum());
573  if (mPointDoubleStats[minPixelZScore] > dValue)
574  mPointDoubleStats[minPixelZScore] = dValue;
575  }
576 
577  Statistics maxPixelZScoreStats = cp->GetStatistic(
578  ControlMeasureLogData::MaximumPixelZScore);
579 
580  if (maxPixelZScoreStats.ValidPixels()) {
581  dValue = fabs(maxPixelZScoreStats.Maximum());
582  if (mPointDoubleStats[maxPixelZScore] > dValue)
583  mPointDoubleStats[maxPixelZScore] = dValue;
584  }
585  }
586 
587  // Average Residuals
588  mPointDoubleStats[avgResidual] = residualMagStats.Average();
589 
590  // Average Shift
591  mPointDoubleStats[avgPixelShift] = pixelShiftStats.Average();
592  }
593 
594 
595  void ControlNetStatistics::UpdateMinMaxStats(const Statistics & stats,
596  ePointDoubleStats min, ePointDoubleStats max) {
597  if (stats.ValidPixels()) {
598  if (mPointDoubleStats[min] != Null) {
599  mPointDoubleStats[min] = qMin(
600  mPointDoubleStats[min], fabs(stats.Minimum()));
601  }
602  else {
603  mPointDoubleStats[min] = fabs(stats.Minimum());
604  }
605 
606  if (mPointDoubleStats[max] != Null) {
607  mPointDoubleStats[max] = qMax(
608  mPointDoubleStats[max], fabs(stats.Maximum()));
609  }
610  else {
611  mPointDoubleStats[max] = fabs(stats.Maximum());
612  }
613  }
614  }
615 }
616