USGS

Isis 3.0 Object Programmers' Reference

 Home
ProcessPolygons.cpp
1 #include "ProcessPolygons.h"
2 #include "PolygonTools.h"
3 #include "Application.h"
4 #include "BoxcarCachingAlgorithm.h"
5 #include "SpecialPixel.h"
6 
7 #include "geos/geom/CoordinateSequence.h"
8 #include "geos/geom/CoordinateArraySequence.h"
9 #include "geos/geom/Coordinate.h"
10 #include "geos/geom/Envelope.h"
11 #include "geos/geom/LineString.h"
12 #include "geos/geom/Geometry.h"
13 #include "geos/geom/Point.h"
14 #include "geos/geom/prep/PreparedGeometryFactory.h"
15 #include "geos/geom/prep/PreparedPolygon.h"
16 #include "geos/util/IllegalArgumentException.h"
17 #include "geos/operation/overlay/snap/GeometrySnapper.h"
18 
19 using namespace std;
20 namespace Isis {
21 
22  ProcessPolygons::ProcessPolygons() {
23  m_useCenter = true;
24  m_average = NULL;
25  m_count = NULL;
26  m_imagePoly = NULL;
27  }
28 
29 
30 
38  void ProcessPolygons::Rasterize(std::vector<double> &samples,
39  std::vector<double> &lines,
40  std::vector<double> &values) {
41 
42  m_sampleVertices = samples;
43  m_lineVertices = lines;
44  m_dns = values;
45  FillPolygon(0);
46  }
47 
48 
49 
59  void ProcessPolygons::Rasterize(std::vector<double> &samples,
60  std::vector<double> &lines,
61  int &band, double &value) {
62 
63  m_sampleVertices = samples;
64  m_lineVertices = lines;
65  m_band = band;
66  m_dns.clear();
67  m_dns.push_back(value);
68 
69  FillPolygon(1);
70  }
71 
72 
83  void ProcessPolygons::FillPolygon(int Flag) {
84 
85  // Create a sample/line polygon for the input pixel vertices
86  geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateArraySequence();
87  for (unsigned int i = 0; i < m_sampleVertices.size(); i++) {
88  pts->add(geos::geom::Coordinate(m_sampleVertices[i], m_lineVertices[i]));
89  }
90  pts->add(geos::geom::Coordinate(m_sampleVertices[0], m_lineVertices[0]));
91 
92  try {
93  // Create a polygon from the pixel vertices. This polygon may have spikes or other
94  // problems such as multiple polygons. Despike, then make sure we have a single polygon.
95  // Do not rasterize pixel if despiking fails or there are multiple polygons.
96  geos::geom::Polygon *spikedPixelPoly = Isis::globalFactory.createPolygon(
97  globalFactory.createLinearRing(pts), NULL);
98 
99  const geos::geom::Polygon *projectedInputPixelPoly;
100 
101  if (spikedPixelPoly->isValid()) {
102  projectedInputPixelPoly = spikedPixelPoly;
103  }
104  else {
105  geos::geom::MultiPolygon *despikedPixelPoly;
106  try {
107  despikedPixelPoly = PolygonTools::Despike(spikedPixelPoly);
108  }
109  catch (IException &e) {
110  delete spikedPixelPoly;
111  return;
112  }
113 
114  try {
115  despikedPixelPoly = PolygonTools::Despike(spikedPixelPoly);
116  }
117  catch (IException &e) {
118  delete spikedPixelPoly;
119  return;
120  }
121 
122  if (despikedPixelPoly->getNumGeometries() > 1) return;
123  projectedInputPixelPoly =
124  dynamic_cast<const geos::geom::Polygon *>(despikedPixelPoly->getGeometryN(0));
125  }
126 
127  /* If there is not an intersecting polygon, there is no reason to go on.*/
128  if (!projectedInputPixelPoly->intersects(m_imagePoly)) return;
129 
130  geos::geom::MultiPolygon *intersectPoly = PolygonTools::MakeMultiPolygon(
131  m_imagePoly->intersection(projectedInputPixelPoly));
132  geos::geom::prep::PreparedPolygon *preparedPoly =
133  new geos::geom::prep::PreparedPolygon(intersectPoly);
134  const geos::geom::Envelope *envelope = intersectPoly->getEnvelopeInternal();
135 
136  // Go thru each coord. in the envelope and ask if it is within the polygon
137  for (double x = floor(envelope->getMinX()); x <= ceil(envelope->getMaxX()); x++) {
138  if (x == 0) continue;
139 
140  for (double y = floor(envelope->getMinY()); y <= ceil(envelope->getMaxY()); y++) {
141  if (y == 0) continue;
142 
143  bool contains;
144 
145  if (m_useCenter) {
146  geos::geom::Coordinate c(x, y);
147  geos::geom::Point *p = Isis::globalFactory.createPoint(c);
148  contains = preparedPoly->contains(p);
149  delete p;
150  }
151  else {
152  geos::geom::CoordinateSequence *tpts = new geos::geom::CoordinateArraySequence();
153  tpts->add(geos::geom::Coordinate(x - 0.5, y - 0.5));
154  tpts->add(geos::geom::Coordinate(x + 0.5, y - 0.5));
155  tpts->add(geos::geom::Coordinate(x + 0.5, y + 0.5));
156  tpts->add(geos::geom::Coordinate(x - 0.5, y + 0.5));
157  tpts->add(geos::geom::Coordinate(x - 0.5, y - 0.5));
158 
159  geos::geom::Polygon *outPixelFootPrint = Isis::globalFactory.createPolygon(
160  globalFactory.createLinearRing(tpts), NULL);
161  contains = preparedPoly->intersects(outPixelFootPrint);
162  delete outPixelFootPrint;
163  }
164 
165  if (contains) {
166 
167  // Read spectral noodle from samp, line position
168  m_average->SetBasePosition((int)(x + 0.5), (int)(y + 0.5), 1);
169  this->OutputCubes[0]->read(*m_average);
170  m_count->SetBasePosition((int)(x + 0.5), (int)(y + 0.5), 1);
171  this->OutputCubes[1]->read(*m_count);
172 
173  // Process each band in the buffer
174  for (unsigned int i = 0; i < m_dns.size(); i++) {
175 
176  int band;
177  if (Flag == 0) {
178  band = i;
179  }
180  else {
181  band = m_band - 1;
182  }
183 
184  double inputDn = m_dns[i];
185 
186  // The input dn is good
187  if (IsValidPixel(inputDn)) {
188  if (IsValidPixel((*m_average)[band])) {
189  double currentCount = (*m_count)[band];
190  double newCount = ++((*m_count)[band]);
191  double currentAverage = (*m_average)[band];
192  (*m_average)[band] = ((currentAverage * currentCount) + inputDn) / newCount;
193  }
194  else {
195  (*m_average)[band] = inputDn;
196  (*m_count)[band] = 1;
197  }
198  }
199 
200  // The input dn is special
201  else {
202  if (((*m_average)[band] == Isis::Null) || (inputDn != Null)) {
203  (*m_average)[band] = inputDn;
204  }
205  }
206 
207  } /*End for each band*/
208 
209  // Write spectral noodles back out to average and count cubes
210  this->OutputCubes[0]->write(*m_average);
211  this->OutputCubes[1]->write(*m_count);
212 
213  } /*End if (contains)*/
214 
215  } /*End for y*/
216 
217  } /*End for x*/
218 
219  delete projectedInputPixelPoly;
220  delete intersectPoly;
221  delete preparedPoly;
222 
223  } /*end try*/
224 
225  catch(geos::util::IllegalArgumentException *ill) {
226  QString msg = "ERROR! geos exception 1 [";
227  msg += (QString)ill->what() + "]";
228  delete ill;
229  throw IException(IException::Programmer, msg, _FILEINFO_);
230  }/*end catch*/
231 
232  }
233 
234 
239  void ProcessPolygons::EndProcess() {
240 
241  delete m_imagePoly;
242  delete m_average;
243  delete m_count;
244  Process::EndProcess();
245  }
246 
247 
252  void ProcessPolygons::Finalize() {
253 
254  delete m_imagePoly;
255  delete m_average;
256  delete m_count;
257  Process::Finalize();
258  }
259 
260 
269  Isis::Cube *ProcessPolygons::AppendOutputCube(const QString &avgFileName,
270  const QString &countFileName) {
271 
272  FileName *file = new FileName(avgFileName);
273  QString path = file->path();
274  QString filename = file->baseName();
275  QString extension = file->extension();
276 
277  /*Open the average file with read/write permission*/
278  Cube *averageCube = new Cube();
279  averageCube->open(avgFileName, "rw");
280  AddOutputCube(averageCube);
281 
282  /*Now open the count file with read/write permission*/
283  Cube *countCube = new Cube();
284 
285  if (countFileName == "") {
286  /*if the countFileName was set to nothing, then we use the default count
287  file name.*/
288  QString openFile = path + "/" + filename + "-count-." + extension;
289  countCube->open(openFile, "rw");
290 
291  }
292  else {
293  countCube->open(countFileName, "rw");
294  }
295 
296  AddOutputCube(countCube);
297  return countCube;
298  }
299 
300 
301 
311  void ProcessPolygons::SetStatCubes(const QString &avgFileName, const
312  QString &countFileName,
313  Isis::CubeAttributeOutput &atts,
314  const int nsamps, const int nlines,
315  const int nbands) {
316 
317  this->Process::SetOutputCube(avgFileName, atts, nsamps, nlines, nbands);
318  this->Process::SetOutputCube(countFileName, atts, nsamps, nlines, nbands);
319 
320  OutputCubes[0]->addCachingAlgorithm(new BoxcarCachingAlgorithm());
321  OutputCubes[1]->addCachingAlgorithm(new BoxcarCachingAlgorithm());
322 
323  geos::geom::CoordinateArraySequence imagePts;
324 
325  imagePts.add(geos::geom::Coordinate(0.0, 0.0));
326  imagePts.add(geos::geom::Coordinate(0.0, this->OutputCubes[0]->lineCount()));
327  imagePts.add(geos::geom::Coordinate(this->OutputCubes[0]->sampleCount(),
328  this->OutputCubes[0]->lineCount()));
329  imagePts.add(geos::geom::Coordinate(this->OutputCubes[0]->sampleCount(), 0.0));
330  imagePts.add(geos::geom::Coordinate(0.0, 0.0));
331 
332  m_imagePoly = Isis::globalFactory.createPolygon(
333  globalFactory.createLinearRing(imagePts), NULL);
334 
335  m_average = new Brick(*this->OutputCubes[0], 1, 1, nbands);
336  m_count = new Brick(*this->OutputCubes[1], 1, 1, nbands);
337  }
338 
339 
340 
349  void ProcessPolygons::SetStatCubes(const QString &parameter,
350  const int nsamps, const int nlines,
351  const int nbands) {
352 
353  QString avgString =
354  Application::GetUserInterface().GetFileName(parameter);
355 
356  Isis::CubeAttributeOutput atts =
357  Application::GetUserInterface().GetOutputAttribute(parameter);
358 
359  FileName *file = new FileName(avgString);
360  QString path = file->path();
361  QString filename = file->baseName();
362  QString countString = path + "/" + filename + "-count";
363  SetStatCubes(avgString, countString, atts, nsamps, nlines, nbands);
364 
365  }
366 
367 
374  void ProcessPolygons::SetIntersectAlgorithm(const bool useCenter) {
375  m_useCenter = useCenter;
376  }
377 
378 
379 } /* end namespace isis*/
380