Isis 3 Programmer Reference
GridPolygonSeeder.cpp
Go to the documentation of this file.
1 
24 #include <string>
25 #include <vector>
26 #include <cmath>
27 
28 #include "Pvl.h"
29 #include "PvlGroup.h"
30 #include "IException.h"
31 #include "PolygonTools.h"
32 #include "IString.h"
33 
34 #include "GridPolygonSeeder.h"
35 
36 namespace Isis {
37 
46  Parse(pvl);
47  };
48 
49 
68  std::vector<geos::geom::Point *> GridPolygonSeeder::Seed(const geos::geom::MultiPolygon *lonLatPoly) {
69  //Projection *proj) {
70  /*if (proj == NULL) {
71  QString msg = "No Projection object available";
72  throw iException::Message(iException::Programmer, msg, _FILEINFO_);
73  }*/
74 
75  if(!p_subGrid)
76  //return SeedGrid(lonLatPoly, proj);
77  return SeedGrid(lonLatPoly);
78  else
79  //return SeedSubGrid(lonLatPoly, proj);
80  return SeedSubGrid(lonLatPoly);
81 
82  }
83 
84  std::vector<geos::geom::Point *> GridPolygonSeeder::SeedGrid(const geos::geom::MultiPolygon *multiPoly) {
85 
86  // Storage for the points to be returned
87  std::vector<geos::geom::Point *> points;
88 
89  // Create some things we will need shortly
90  const geos::geom::Envelope *polyBoundBox = multiPoly->getEnvelopeInternal();
91 
92  // Call the parents standardTests member
93  QString msg = StandardTests(multiPoly, polyBoundBox);
94  if(!msg.isEmpty()) {
95  return points;
96  }
97 
98  // Do grid specific tests to make sure this poly should be seeded
99  // (none for now)
100 
101  // Starting at the centroid of the xy polygon populate the polygon with a
102  // grid of points with the requested spacing
103  geos::geom::Point *centroid = multiPoly->getCentroid();
104  double centerX = centroid->getX();
105  double centerY = centroid->getY();
106  delete centroid;
107 
108  int xStepsLeft = (int)((centerX - polyBoundBox->getMinX()) / p_Xspacing + 0.5);
109  int yStepsLeft = (int)((centerY - polyBoundBox->getMinY()) / p_Yspacing + 0.5);
110  double dRealMinX = centerX - (xStepsLeft * p_Xspacing);
111  double dRealMinY = centerY - (yStepsLeft * p_Yspacing);
112 
113  for(double y = dRealMinY; y <= polyBoundBox->getMaxY(); y += p_Yspacing) {
114  for(double x = dRealMinX; x <= polyBoundBox->getMaxX(); x += p_Xspacing) {
115  geos::geom::Coordinate c(x, y);
116  geos::geom::Point *p = Isis::globalFactory.createPoint(c);
117 
118  if(p->within(multiPoly)) {
119  points.push_back(Isis::globalFactory.createPoint(c));
120  }
121  else {
122  delete p;
123  }
124  }
125  }
126 
127  return points;
128  }
129 
140  std::vector<geos::geom::Point *> GridPolygonSeeder::SeedSubGrid(const geos::geom::MultiPolygon *multiPoly) {
141  //Projection *proj) {
142  // Storage for the points to be returned
143  std::vector<geos::geom::Point *> points;
144 
145  // Create some things we will need shortly
146  //geos::geom::MultiPolygon *xymp = PolygonTools::LatLonToXY(*lonLatPoly, proj);
147  const geos::geom::Envelope *polyBoundBox = multiPoly->getEnvelopeInternal();
148 
149  // Call the parents standardTests member
150  QString msg = StandardTests(multiPoly, polyBoundBox);
151  if(!msg.isEmpty()) {
152  return points;
153  }
154 
155  geos::geom::Point *centroid = multiPoly->getCentroid();
156  double centerX = centroid->getX();
157  double centerY = centroid->getY();
158  delete centroid;
159 
160  // Do grid specific tests to make sure this poly should be seeded
161  // (none for now)
162 
171  enum PointStatus {
172  pointShouldCheck,
173  pointShouldSubGridCheck,
174  pointFound,
175  pointNotFound,
176  pointCantFind
177  };
178 
179  // For maintaining an idea of what's going on in this polygon, we needs to know the dimensions
180  // of the grid.
181  int xSteps = (int)((polyBoundBox->getMaxX() - polyBoundBox->getMinX()) / p_Xspacing + 1.5);
182  int ySteps = (int)((polyBoundBox->getMaxY() - polyBoundBox->getMinY()) / p_Yspacing + 1.5);
183  PointStatus pointCheck[xSteps][ySteps];
184 
185  // Initialize our grid of point status'
186  for(int y = 0; y < ySteps; y++) {
187  for(int x = 0; x < xSteps; x++) {
188  pointCheck[x][y] = pointShouldCheck;
189  }
190  }
191 
200  int precision = (int)pow(0.5 / MinimumThickness(), 0.5) * 2;
201  bool bGridCleared = true;
202  int xStepsToCentroid = (int)((centerX - polyBoundBox->getMinX()) / p_Xspacing + 0.5);
203  int yStepsToCentroid = (int)((centerY - polyBoundBox->getMinY()) / p_Yspacing + 0.5);
204  double dRealMinX = centerX - (xStepsToCentroid * p_Xspacing);
205  double dRealMinY = centerY - (yStepsToCentroid * p_Yspacing);
206 
207  do {
208  // gridCleared is true if we did nothing, if we performed any actions on the grid
209  // it becomes false and another pass should be used.
210  bGridCleared = true;
211 
212  for(int y = 0; y < ySteps; y++) {
213  double centerY = dRealMinY + p_Yspacing * y;
214  for(int x = 0; x < xSteps; x++) {
215  double centerX = dRealMinX + p_Xspacing * x;
216  geos::geom::Point *p = NULL;
217 
218  // pointShouldCheck tells us we need to check center. Calling
219  // CheckSubGrid with precision=0 will do this for us.
220  if(pointCheck[x][y] == pointShouldCheck) {
221  p = CheckSubGrid(*multiPoly, centerX, centerY, 0);
222  }
223  // pointShouldSubGridCheck tells us we're next to a grid
224  // square where a point was found, so check in depth
225  else if(pointCheck[x][y] == pointShouldSubGridCheck) {
226  p = CheckSubGrid(*multiPoly, centerX, centerY, precision);
227  }
228 
229  // If we found a point, verify we can setCoordinate and save the point
230  if(p != NULL) {
231  // Convert the x/y point to a lon/lat point
232  /*if (proj->SetCoordinate(p->getX(),p->getY())) {
233  points.push_back(Isis::globalFactory.createPoint(
234  geos::geom::Coordinate(proj->UniversalLongitude(),
235  proj->UniversalLatitude())));
236  }
237  else {
238  IString msg = "Unable to convert [(" + IString(x) + ",";
239  msg += IString(y) + ")] to a (lon,lat)";
240  throw iException::Message(iException::Programmer, msg, _FILEINFO_);
241  }*/
242  points.push_back(Isis::globalFactory.createPoint(
243  geos::geom::Coordinate(p->getX(), p->getY())));
244 
245  // We found something new and need a new pass
246  bGridCleared = false;
247  pointCheck[x][y] = pointFound;
248  }
249  else {
250  if(pointCheck[x][y] == pointShouldCheck) {
251  pointCheck[x][y] = pointNotFound;
252  }
253  else if(pointCheck[x][y] == pointShouldSubGridCheck) {
254  pointCheck[x][y] = pointCantFind;
255  }
256  }
257  }
258  }
259 
260  // now that the grid has been updated with it's founds, we can look for subgrid checks
261  for(int y = 0; y < ySteps; y++) {
262  for(int x = 0; x < xSteps; x++) {
263  if(pointCheck[x][y] == pointFound) {
264  for(int yOff = -1; yOff <= 1; yOff++) {
265  for(int xOff = -1; xOff <= 1; xOff++) {
266  if(x + xOff >= 0 && x + xOff < xSteps &&
267  y + yOff >= 0 && y + yOff < ySteps &&
268  pointCheck[x+xOff][y+yOff] == pointNotFound) {
269 
270  pointCheck[x+xOff][y+yOff] = pointShouldSubGridCheck;
271 
272  // We need to do a searchso we need another pass
273  bGridCleared = false;
274  }
275  }
276  }
277  }
278  }
279  }
280 
281  }
282  while(!bGridCleared);
283 
284  return points;
285  }
286 
308  geos::geom::Point *GridPolygonSeeder::CheckSubGrid(const geos::geom::MultiPolygon &xymp, const double &centerX,
309  const double &centerY, const int &precision) {
310  // We'll make a 2D array detailing which points to check, and which not to, in this rectangle.
311  // Figure out how many points across and vertically we need to check
312  int gridSize = 1;
313  for(int prec = 0; prec < precision && prec < 6; prec ++) {
314  // Maybe solve the recurrence relation for a single equation??
315  gridSize = gridSize * 2 + 1;
316  }
317 
318  // These are the possible values in which the 2D array can be. We need the transition value
319  // gridNewCheckPt to not count new points as old points.
320  enum GridPoint {
321  gridEmpty,
322  gridNewCheckPt,
323  gridCheckPt
324  };
325 
326  GridPoint grid[gridSize][gridSize];
327 
328  for(int y = 0; y < gridSize; y++) {
329  for(int x = 0; x < gridSize; x++) {
330  grid[x][y] = gridEmpty;
331  }
332  }
333 
334  // Precision 0: Always center, this is always true
335  grid[gridSize/2][gridSize/2] = gridCheckPt;
336 
337  // now populate the grid with what we wish to check for
338  for(int prec = 0; prec < precision; prec ++) {
339  // This tells us how far over in the 2D array to go at a precision from already found points
340  int checkDist = (gridSize + 1) / (int)(4 * (pow(2.0, prec)) + 0.5);
341 
342  // Search the grid for already found points and set everything checkDist away to be checked too
343  for(int y = 0; y < gridSize; y++) {
344  for(int x = 0; x < gridSize; x++) {
345  if(grid[x][y] == gridCheckPt) {
346  // We should never overwrite found points, the checkDist should assure this wont happen
347  if(x - checkDist > 0) grid[x-checkDist][y] = gridNewCheckPt;
348  if(y - checkDist > 0) grid[x][y-checkDist] = gridNewCheckPt;
349  if(x + checkDist < gridSize) grid[x+checkDist][y] = gridNewCheckPt;
350  if(y + checkDist < gridSize) grid[x][y+checkDist] = gridNewCheckPt;
351  }
352  }
353  }
354 
355  // Convert temporary check pt to final.
356  // Do this separately to avoid changing the data we're processing
357  for(int y = 0; y < gridSize; y++) {
358  for(int x = 0; x < gridSize; x++) {
359  if(grid[x][y] == gridNewCheckPt) grid[x][y] = gridCheckPt;
360  }
361  }
362  }
363 
364  // We now have a grid of points to check inside this grid square. Figure out the
365  // distance each of these subsquare pixels are worth.
366  double deltaXSize = p_Xspacing / (gridSize + 1);
367  double deltaYSize = p_Yspacing / (gridSize + 1);
368 
369  geos::geom::Point *result = NULL;
370  // Now loop through this grid, checking each point that needs checked, return the first valid pt
371  for(int y = 0; !result && y < gridSize; y++) {
372  for(int x = 0; !result && x < gridSize; x++) {
373  if(grid[x][y] != gridCheckPt) continue;
374 
375  double xPos = centerX + (x - gridSize / 2) * deltaXSize;
376  double yPos = centerY + (y - gridSize / 2) * deltaYSize;
377  geos::geom::Coordinate c(xPos, yPos);
378  geos::geom::Point *p = Isis::globalFactory.createPoint(c);
379  if(p->within(&xymp)) {
380  result = p;
381  }
382  else {
383  delete p;
384  }
385  }
386  }
387 
388  return result;
389  }
390 
398  // Call the parents Parse method
400 
401  // Pull parameters specific to this algorithm out
402  try {
403  // Get info from Algorithm group
404  PvlGroup &algo = pvl.findGroup("PolygonSeederAlgorithm", Pvl::Traverse);
405  PvlGroup &invalgo = invalidInput->findGroup("PolygonSeederAlgorithm",
406  Pvl::Traverse);
407 
408  // Set the spacing
409  p_Xspacing = 0.0;
410  if(algo.hasKeyword("XSpacing")) {
411  p_Xspacing = (double) algo["XSpacing"];
412  if(invalgo.hasKeyword("XSpacing")) {
413  invalgo.deleteKeyword("XSpacing");
414  }
415  }
416  else {
417  QString msg = "PVL for GridPolygonSeeder must contain [XSpacing] in [";
418  msg += pvl.fileName() + "]";
420  }
421 
422  p_Yspacing = 0.0;
423  if(algo.hasKeyword("YSpacing")) {
424  p_Yspacing = (double) algo["YSpacing"];
425  if(invalgo.hasKeyword("YSpacing")) {
426  invalgo.deleteKeyword("YSpacing");
427  }
428  }
429  else {
430  QString msg = "PVL for GridPolygonSeeder must contain [YSpacing] in [";
431  msg += pvl.fileName() + "]";
433  }
434 
435  p_subGrid = false;
436  if(algo.hasKeyword("SubGrid")) {
437  p_subGrid = IString((QString)algo["SubGrid"]).UpCase() != "FALSE";
438  if(invalgo.hasKeyword("SubGrid")) {
439  invalgo.deleteKeyword("SubGrid");
440  }
441  }
442  }
443  catch(IException &e) {
444  QString msg = "Improper format for PolygonSeeder PVL [" + pvl.fileName() + "]";
445  throw IException(e, IException::User, msg, _FILEINFO_);
446  }
447 
448  if(p_Xspacing <= 0.0) {
449  IString msg = "X Spacing must be greater that 0.0 [(" + IString(p_Xspacing) + "]";
451  }
452  if(p_Yspacing <= 0.0) {
453  IString msg = "Y Spacing must be greater that 0.0 [(" + IString(p_Yspacing) + "]";
455  }
456  }
457 
459  PvlGroup pluginInfo(grpName);
460 
461  PvlKeyword name("Name", Algorithm());
462  PvlKeyword minThickness("MinimumThickness", toString(MinimumThickness()));
463  PvlKeyword minArea("MinimumArea", toString(MinimumArea()));
464  PvlKeyword xSpac("XSpacing", toString(p_Xspacing));
465  PvlKeyword ySpac("YSpacing", toString(p_Yspacing));
466  PvlKeyword subGrid("SubGrid", toString(p_subGrid));
467 
468  pluginInfo.addKeyword(name);
469  pluginInfo.addKeyword(minThickness);
470  pluginInfo.addKeyword(minArea);
471  pluginInfo.addKeyword(xSpac);
472  pluginInfo.addKeyword(ySpac);
473  pluginInfo.addKeyword(subGrid);
474 
475  return pluginInfo;
476  }
477 
478 }; // End of namespace Isis
479 
480 
493  return new Isis::GridPolygonSeeder(pvl);
494 }
495 
bool hasKeyword(const QString &name) const
Check to see if a keyword exists.
virtual void Parse(Pvl &pvl)
Parse the GridPolygonSeeder spicific parameters from the PVL.
PvlGroupIterator findGroup(const QString &name, PvlGroupIterator beg, PvlGroupIterator end)
Find a group with the specified name, within these indexes.
Definition: PvlObject.h:141
virtual PvlGroup PluginParameters(QString grpName)
Plugin parameters.
QString Algorithm() const
The name of the algorithm, read from the Name Keyword in the PolygonSeeder Pvl passed into the constr...
Pvl * invalidInput
The Pvl passed in by the constructor minus what was used.
Definition: PolygonSeeder.h:94
This class is used as the base class for all PolygonSeeder objects.
Definition: PolygonSeeder.h:63
Search child objects.
Definition: PvlObject.h:170
Seed points using a grid.
QString toString(bool boolToConvert)
Global function to convert a boolean to a string.
Definition: IString.cpp:226
void addKeyword(const PvlKeyword &keyword, const InsertMode mode=Append)
Add a keyword to the container.
Isis::PolygonSeeder * GridPolygonSeederPlugin(Isis::Pvl &pvl)
Create a GridPolygonSeeder object.
Contains multiple PvlContainers.
Definition: PvlGroup.h:57
#define _FILEINFO_
Macro for the filename and line number.
Definition: IException.h:40
A type of error that could only have occurred due to a mistake on the user&#39;s part (e...
Definition: IException.h:142
A single keyword-value pair.
Definition: PvlKeyword.h:98
virtual void Parse(Pvl &pvl)
Initialize parameters in the PolygonSeeder class using a PVL specification.
Container for cube-like labels.
Definition: Pvl.h:135
geos::geom::Point * CheckSubGrid(const geos::geom::MultiPolygon &, const double &, const double &, const int &)
This method is used to search for a valid point, on the polygon, within the square whose center is de...
QString StandardTests(const geos::geom::MultiPolygon *multiPoly, const geos::geom::Envelope *polyBoundBox)
Check the polygon to see if it meets standard criteria.
GridPolygonSeeder(Pvl &pvl)
Construct a GridPolygonSeeder algorithm.
double MinimumThickness()
Return the minimum allowed thickness of the polygon.
QString fileName() const
Returns the filename used to initialise the Pvl object.
Definition: PvlContainer.h:246
double MinimumArea()
Return the minimum allowed area of the polygon.
Isis exception class.
Definition: IException.h:107
Adds specific functionality to C++ strings.
Definition: IString.h:181
Namespace for ISIS/Bullet specific routines.
Definition: Apollo.h:31
std::vector< geos::geom::Point * > Seed(const geos::geom::MultiPolygon *mp)
Seed a polygon with points.
std::vector< geos::geom::Point * > SeedSubGrid(const geos::geom::MultiPolygon *mp)
This method works a lot like SeedGrid, except around the edges of known polygons. ...
IString UpCase()
Converst any lower case characters in the object IString with uppercase characters.
Definition: IString.cpp:632
void deleteKeyword(const QString &name)
Remove a specified keyword.