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  //map< QString, vector<double> >::iterator it;
333  map<QString, bool>::iterator it;
334  // imgSamples, imgLines, imgTotalPoints, imgIgnoredPoints, imgFixedPoints, imgLockedPoints, imgLocked, imgConstrainedPoints, imgFreePoints, imgConvexHullArea, imgConvexHullRatio
335 
336  // Log into the output file
337  ostm << "Filename, SerialNumber, TotalPoints, PointsIgnored, PointsEditLocked, Fixed, Constrained, Free, ConvexHullRatio" << endl;
338  //for (it = mImageMap.begin(); it != mImageMap.end(); it++) {
339  for (it = mSerialNumMap.begin(); it != mSerialNumMap.end(); it++) {
340  ostm << mSerialNumList.FileName((*it).first) << ", " << (*it).first << ", ";
341  bool serialNumExists = (*it).second ;
342  if (serialNumExists) {
343  vector<double>imgStats = mImageMap[(*it).first] ;
344  ostm << imgStats[imgTotalPoints]<< ", " << imgStats[imgIgnoredPoints] << ", " ;
345  ostm << imgStats[imgLockedPoints] << ", " << imgStats[imgFixedPoints] << ", " ;
346  ostm << imgStats[imgConstrainedPoints] << ", " << imgStats[imgFreePoints] << ", ";
347  ostm << imgStats[ imgConvexHullRatio] << endl;
348  }
349  else {
350  ostm << "0, 0, 0, 0, 0, 0, 0" << endl;
351  }
352  }
353  ostm.close();
354  }
355 
356 
366  vector<double> ControlNetStatistics::GetImageStatsBySerialNum(QString psSerialNum) const {
367  return (*mImageMap.find(psSerialNum)).second;
368  }
369 
370 
380  void ControlNetStatistics::GeneratePointStats(const QString &psPointFile) {
381  Isis::FileName outFile(psPointFile);
382 
383  ofstream ostm;
384  QString outName(outFile.expanded());
385  ostm.open(outName.toAscii().data(), std::ios::out);
386  ostm << " PointId, PointType, PointIgnore, PointEditLock, TotalMeasures, MeasuresValid, MeasuresIgnore, MeasuresEditLock," << endl;
387 
388  int iNumPoints = mCNet->GetNumPoints();
389 
390  // Initialise the Progress object
391  if (mProgress != NULL && iNumPoints > 0) {
392  mProgress->SetText("Point Stats: Loading Control Points...");
393  mProgress->SetMaximumSteps(iNumPoints);
394  mProgress->CheckStatus();
395  }
396 
397  for (int i = 0; i < iNumPoints; i++) {
398  const ControlPoint *cPoint = mCNet->GetPoint(i);
399  int iNumMeasures = cPoint->GetNumMeasures();
400  int iValidMeasures = cPoint->GetNumValidMeasures();
401  int iIgnoredMeasures = iNumMeasures - iValidMeasures;
402 
403  // Log into the output file
404  ostm << cPoint->GetId() << ", " << sPointType[(int)cPoint->GetType()] << ", " << sBoolean[(int)cPoint->IsIgnored()] << ", " ;
405  ostm << sBoolean[(int)cPoint->IsEditLocked()] << ", " << iNumMeasures << ", " << iValidMeasures << ", ";
406  ostm << iIgnoredMeasures << ", " << cPoint->GetNumLockedMeasures() << endl;
407 
408  // Update Progress
409  if (mProgress != NULL)
410  mProgress->CheckStatus();
411  }
412  ostm.close();
413  }
414 
415 
421  void ControlNetStatistics::GetPointIntStats() {
422  // Init all the entries
423  // totalPoints, validPoints, ignoredPoints, fixedPoints, constrainedPoints, editLockedPoints,
424  // totalMeasures, validMeasures, ignoredMeasures, editLockedMeasures
425  for (int i=0; i<numPointIntStats; i++) {
426  mPointIntStats[i] = 0;
427  }
428 
429  int iNumPoints = mCNet->GetNumPoints();
430 
431  // totalPoints
432  mPointIntStats[totalPoints] = iNumPoints;
433 
434  for (int i = 0; i < iNumPoints; i++) {
435  if (!mCNet->GetPoint(i)->IsIgnored()) {
436  // validPoints
437  mPointIntStats[validPoints]++;
438  }
439  else {
440  // ignoredPoints
441  mPointIntStats[ignoredPoints]++;
442  }
443 
444  // fixedPoints
445  if (mCNet->GetPoint(i)->GetType() == ControlPoint::Fixed)
446  mPointIntStats[fixedPoints]++;
447 
448  // constrainedPoints
449  if (mCNet->GetPoint(i)->GetType() == ControlPoint::Constrained)
450  mPointIntStats[constrainedPoints]++;
451 
452  // free points
453  if (mCNet->GetPoint(i)->GetType() == ControlPoint::Free)
454  mPointIntStats[freePoints]++;
455 
456  // editLockedPoints
457  if (mCNet->GetPoint(i)->IsEditLocked()) {
458  mPointIntStats[editLockedPoints]++;
459  }
460 
461  // totalMeasures
462  mPointIntStats[totalMeasures] += mCNet->GetPoint(i)->GetNumMeasures();
463 
464  // validMeasures
465  mPointIntStats[validMeasures] += mCNet->GetPoint(i)->GetNumValidMeasures();
466 
467  // editLockedMeasures
468  mPointIntStats[editLockedMeasures] += mCNet->GetPoint(i)->GetNumLockedMeasures();
469  }
470  // ignoredMeasures
471  mPointIntStats[ignoredMeasures] = mPointIntStats[totalMeasures] - mPointIntStats[validMeasures];
472  }
473 
474 
480  void ControlNetStatistics::InitPointDoubleStats() {
481  for (int i = 0; i < numPointDblStats; i++) {
482  mPointDoubleStats[i] = Null;
483  }
484  }
485 
486 
493  void ControlNetStatistics::GetPointDoubleStats() {
494  InitPointDoubleStats();
495 
496  int iNumPoints = mCNet->GetNumPoints();
497  double dValue = 0;
498 
499  Statistics residualMagStats;
500  Statistics pixelShiftStats;
501 
502  for (int i = 0; i < iNumPoints; i++) {
503  ControlPoint * cp = mCNet->GetPoint(i);
504 
505  if (!cp->IsIgnored()) {
506  for (int cmIndex = 0; cmIndex < cp->GetNumMeasures(); cmIndex++) {
507  ControlMeasure *cm = cp->GetMeasure(cmIndex);
508 
509  if (!cm->IsIgnored()) {
510  residualMagStats.AddData(cm->GetResidualMagnitude());
511 
512  if (!IsSpecial(cm->GetPixelShift()))
513  pixelShiftStats.AddData(fabs(cm->GetPixelShift()));
514  }
515  }
516  }
517 
518  Statistics resMagStats = cp->GetStatistic(
519  &ControlMeasure::GetResidualMagnitude);
520  UpdateMinMaxStats(resMagStats, minResidual, maxResidual);
521 
522  Statistics resLineStats = cp->GetStatistic(
523  &ControlMeasure::GetLineResidual);
524  UpdateMinMaxStats(resLineStats, minLineResidual, maxLineResidual);
525 
526  Statistics resSampStats = cp->GetStatistic(
527  &ControlMeasure::GetSampleResidual);
528  UpdateMinMaxStats(resSampStats, minSampleResidual, maxSampleResidual);
529 
530  Statistics pixShiftStats = cp->GetStatistic(
531  &ControlMeasure::GetPixelShift);
532  UpdateMinMaxStats(pixShiftStats, minPixelShift, maxPixelShift);
533 
534  Statistics lineShiftStats = cp->GetStatistic(
535  &ControlMeasure::GetLineShift);
536  UpdateMinMaxStats(lineShiftStats, minLineShift, maxLineShift);
537 
538  Statistics sampShiftStats = cp->GetStatistic(
539  &ControlMeasure::GetSampleShift);
540  UpdateMinMaxStats(sampShiftStats, minSampleShift, maxSampleShift);
541 
542  Statistics gFitStats = cp->GetStatistic(
543  ControlMeasureLogData::GoodnessOfFit);
544  UpdateMinMaxStats(gFitStats, minGFit, maxGFit);
545 
546  Statistics minPixelZScoreStats = cp->GetStatistic(
547  ControlMeasureLogData::MinimumPixelZScore);
548 
549  if (minPixelZScoreStats.ValidPixels()) {
550  dValue = fabs(minPixelZScoreStats.Minimum());
551  if (mPointDoubleStats[minPixelZScore] > dValue)
552  mPointDoubleStats[minPixelZScore] = dValue;
553  }
554 
555  Statistics maxPixelZScoreStats = cp->GetStatistic(
556  ControlMeasureLogData::MaximumPixelZScore);
557 
558  if (maxPixelZScoreStats.ValidPixels()) {
559  dValue = fabs(maxPixelZScoreStats.Maximum());
560  if (mPointDoubleStats[maxPixelZScore] > dValue)
561  mPointDoubleStats[maxPixelZScore] = dValue;
562  }
563  }
564 
565  // Average Residuals
566  mPointDoubleStats[avgResidual] = residualMagStats.Average();
567 
568  // Average Shift
569  mPointDoubleStats[avgPixelShift] = pixelShiftStats.Average();
570  }
571 
572 
573  void ControlNetStatistics::UpdateMinMaxStats(const Statistics & stats,
574  ePointDoubleStats min, ePointDoubleStats max) {
575  if (stats.ValidPixels()) {
576  if (mPointDoubleStats[min] != Null) {
577  mPointDoubleStats[min] = qMin(
578  mPointDoubleStats[min], fabs(stats.Minimum()));
579  }
580  else {
581  mPointDoubleStats[min] = fabs(stats.Minimum());
582  }
583 
584  if (mPointDoubleStats[max] != Null) {
585  mPointDoubleStats[max] = qMax(
586  mPointDoubleStats[max], fabs(stats.Maximum()));
587  }
588  else {
589  mPointDoubleStats[max] = fabs(stats.Maximum());
590  }
591  }
592  }
593 }
594