Isis 3 Programmer Reference
ProcessPolygons.cpp
1 
6 /* SPDX-License-Identifier: CC0-1.0 */
7 #include "ProcessPolygons.h"
8 #include "PolygonTools.h"
9 #include "Application.h"
10 #include "BoxcarCachingAlgorithm.h"
11 #include "SpecialPixel.h"
12 
13 #include "geos/geom/CoordinateSequence.h"
14 #include "geos/geom/CoordinateArraySequence.h"
15 #include "geos/geom/Coordinate.h"
16 #include "geos/geom/Envelope.h"
17 #include "geos/geom/LineString.h"
18 #include "geos/geom/Geometry.h"
19 #include "geos/geom/Point.h"
20 #include "geos/geom/prep/PreparedGeometryFactory.h"
21 #include "geos/geom/prep/PreparedPolygon.h"
22 #include "geos/util/IllegalArgumentException.h"
23 #include "geos/operation/overlay/snap/GeometrySnapper.h"
24 
25 using namespace std;
26 namespace Isis {
27 
28  ProcessPolygons::ProcessPolygons() {
29  m_useCenter = true;
30  m_average = NULL;
31  m_count = NULL;
32  m_imagePoly = NULL;
33  }
34 
35 
36 
44  void ProcessPolygons::Rasterize(std::vector<double> &samples,
45  std::vector<double> &lines,
46  std::vector<double> &values) {
47 
48  m_sampleVertices = samples;
49  m_lineVertices = lines;
50  m_dns = values;
51  FillPolygon(0);
52  }
53 
54 
55 
65  void ProcessPolygons::Rasterize(std::vector<double> &samples,
66  std::vector<double> &lines,
67  int &band, double &value) {
68 
69  m_sampleVertices = samples;
70  m_lineVertices = lines;
71  m_band = band;
72  m_dns.clear();
73  m_dns.push_back(value);
74 
75  FillPolygon(1);
76  }
77 
78 
89  void ProcessPolygons::FillPolygon(int Flag) {
90 
91  // Create a sample/line polygon for the input pixel vertices
92  geos::geom::CoordinateSequence *pts = new geos::geom::CoordinateArraySequence();
93  for (unsigned int i = 0; i < m_sampleVertices.size(); i++) {
94  pts->add(geos::geom::Coordinate(m_sampleVertices[i], m_lineVertices[i]));
95  }
96  pts->add(geos::geom::Coordinate(m_sampleVertices[0], m_lineVertices[0]));
97 
98  try {
99  // Create a polygon from the pixel vertices. This polygon may have spikes or other
100  // problems such as multiple polygons. Despike, then make sure we have a single polygon.
101  // Do not rasterize pixel if despiking fails or there are multiple polygons.
102  geos::geom::Polygon *spikedPixelPoly = Isis::globalFactory->createPolygon(
103  globalFactory->createLinearRing(pts), NULL);
104 
105  const geos::geom::Polygon *projectedInputPixelPoly;
106 
107  if (spikedPixelPoly->isValid()) {
108  projectedInputPixelPoly = spikedPixelPoly;
109  }
110  else {
111  geos::geom::MultiPolygon *despikedPixelPoly;
112  try {
113  despikedPixelPoly = PolygonTools::Despike(spikedPixelPoly);
114  }
115  catch (IException &e) {
116  delete spikedPixelPoly;
117  return;
118  }
119 
120  try {
121  despikedPixelPoly = PolygonTools::Despike(spikedPixelPoly);
122  }
123  catch (IException &e) {
124  delete spikedPixelPoly;
125  return;
126  }
127 
128  if (despikedPixelPoly->getNumGeometries() > 1) return;
129  projectedInputPixelPoly =
130  dynamic_cast<const geos::geom::Polygon *>(despikedPixelPoly->getGeometryN(0));
131  }
132 
133  /* If there is not an intersecting polygon, there is no reason to go on.*/
134  if (!projectedInputPixelPoly->intersects(m_imagePoly)) return;
135 
136  geos::geom::MultiPolygon *intersectPoly = PolygonTools::MakeMultiPolygon(
137  m_imagePoly->intersection(projectedInputPixelPoly));
138  geos::geom::prep::PreparedPolygon *preparedPoly =
139  new geos::geom::prep::PreparedPolygon(intersectPoly);
140  const geos::geom::Envelope *envelope = intersectPoly->getEnvelopeInternal();
141 
142  // Go thru each coord. in the envelope and ask if it is within the polygon
143  for (double x = floor(envelope->getMinX()); x <= ceil(envelope->getMaxX()); x++) {
144  if (x == 0) continue;
145 
146  for (double y = floor(envelope->getMinY()); y <= ceil(envelope->getMaxY()); y++) {
147  if (y == 0) continue;
148 
149  bool contains;
150 
151  if (m_useCenter) {
152  geos::geom::Coordinate c(x, y);
153  geos::geom::Point *p = Isis::globalFactory->createPoint(c);
154  contains = preparedPoly->contains(p);
155  delete p;
156  }
157  else {
158  geos::geom::CoordinateSequence *tpts = new geos::geom::CoordinateArraySequence();
159  tpts->add(geos::geom::Coordinate(x - 0.5, y - 0.5));
160  tpts->add(geos::geom::Coordinate(x + 0.5, y - 0.5));
161  tpts->add(geos::geom::Coordinate(x + 0.5, y + 0.5));
162  tpts->add(geos::geom::Coordinate(x - 0.5, y + 0.5));
163  tpts->add(geos::geom::Coordinate(x - 0.5, y - 0.5));
164 
165  geos::geom::Polygon *outPixelFootPrint = Isis::globalFactory->createPolygon(
166  globalFactory->createLinearRing(tpts), NULL);
167  contains = preparedPoly->intersects(outPixelFootPrint);
168  delete outPixelFootPrint;
169  }
170 
171  if (contains) {
172 
173  // Read spectral noodle from samp, line position
174  m_average->SetBasePosition((int)(x + 0.5), (int)(y + 0.5), 1);
175  this->OutputCubes[0]->read(*m_average);
176  m_count->SetBasePosition((int)(x + 0.5), (int)(y + 0.5), 1);
177  this->OutputCubes[1]->read(*m_count);
178 
179  // Process each band in the buffer
180  for (unsigned int i = 0; i < m_dns.size(); i++) {
181 
182  int band;
183  if (Flag == 0) {
184  band = i;
185  }
186  else {
187  band = m_band - 1;
188  }
189 
190  double inputDn = m_dns[i];
191 
192  // The input dn is good
193  if (IsValidPixel(inputDn)) {
194  if (IsValidPixel((*m_average)[band])) {
195  double currentCount = (*m_count)[band];
196  double newCount = ++((*m_count)[band]);
197  double currentAverage = (*m_average)[band];
198  (*m_average)[band] = ((currentAverage * currentCount) + inputDn) / newCount;
199  }
200  else {
201  (*m_average)[band] = inputDn;
202  (*m_count)[band] = 1;
203  }
204  }
205 
206  // The input dn is special
207  else {
208  if (((*m_average)[band] == Isis::Null) || (inputDn != Null)) {
209  (*m_average)[band] = inputDn;
210  }
211  }
212 
213  } /*End for each band*/
214 
215  // Write spectral noodles back out to average and count cubes
216  this->OutputCubes[0]->write(*m_average);
217  this->OutputCubes[1]->write(*m_count);
218 
219  } /*End if (contains)*/
220 
221  } /*End for y*/
222 
223  } /*End for x*/
224 
225  delete projectedInputPixelPoly;
226  delete intersectPoly;
227  delete preparedPoly;
228 
229  } /*end try*/
230 
231  catch(geos::util::IllegalArgumentException *ill) {
232  QString msg = "ERROR! geos exception 1 [";
233  msg += (QString)ill->what() + "]";
234  delete ill;
235  throw IException(IException::Programmer, msg, _FILEINFO_);
236  }/*end catch*/
237 
238  }
239 
240 
245  void ProcessPolygons::EndProcess() {
246 
247  delete m_imagePoly;
248  delete m_average;
249  delete m_count;
250  Process::EndProcess();
251  }
252 
253 
258  void ProcessPolygons::Finalize() {
259 
260  delete m_imagePoly;
261  delete m_average;
262  delete m_count;
263  Process::Finalize();
264  }
265 
266 
275  Isis::Cube *ProcessPolygons::AppendOutputCube(const QString &avgFileName,
276  const QString &countFileName) {
277 
278  FileName *file = new FileName(avgFileName);
279  QString path = file->path();
280  QString filename = file->baseName();
281  QString extension = file->extension();
282 
283  /*Open the average file with read/write permission*/
284  Cube *averageCube = new Cube();
285  averageCube->open(avgFileName, "rw");
286  AddOutputCube(averageCube);
287 
288  /*Now open the count file with read/write permission*/
289  Cube *countCube = new Cube();
290 
291  if (countFileName == "") {
292  /*if the countFileName was set to nothing, then we use the default count
293  file name.*/
294  QString openFile = path + "/" + filename + "-count-." + extension;
295  countCube->open(openFile, "rw");
296 
297  }
298  else {
299  countCube->open(countFileName, "rw");
300  }
301 
302  AddOutputCube(countCube);
303  return countCube;
304  }
305 
306 
307 
317  void ProcessPolygons::SetStatCubes(const QString &avgFileName, const
318  QString &countFileName,
319  Isis::CubeAttributeOutput &atts,
320  const int nsamps, const int nlines,
321  const int nbands) {
322 
323  this->Process::SetOutputCube(avgFileName, atts, nsamps, nlines, nbands);
324  this->Process::SetOutputCube(countFileName, atts, nsamps, nlines, nbands);
325 
326  OutputCubes[0]->addCachingAlgorithm(new BoxcarCachingAlgorithm());
327  OutputCubes[1]->addCachingAlgorithm(new BoxcarCachingAlgorithm());
328 
329  geos::geom::CoordinateArraySequence imagePts;
330 
331  imagePts.add(geos::geom::Coordinate(0.0, 0.0));
332  imagePts.add(geos::geom::Coordinate(0.0, this->OutputCubes[0]->lineCount()));
333  imagePts.add(geos::geom::Coordinate(this->OutputCubes[0]->sampleCount(),
334  this->OutputCubes[0]->lineCount()));
335  imagePts.add(geos::geom::Coordinate(this->OutputCubes[0]->sampleCount(), 0.0));
336  imagePts.add(geos::geom::Coordinate(0.0, 0.0));
337 
338  m_imagePoly = Isis::globalFactory->createPolygon(
339  globalFactory->createLinearRing(imagePts), NULL);
340 
341  m_average = new Brick(*this->OutputCubes[0], 1, 1, nbands);
342  m_count = new Brick(*this->OutputCubes[1], 1, 1, nbands);
343  }
344 
345 
346 
355  void ProcessPolygons::SetStatCubes(const QString &parameter,
356  const int nsamps, const int nlines,
357  const int nbands) {
358 
359  QString avgString =
360  Application::GetUserInterface().GetFileName(parameter);
361 
362  Isis::CubeAttributeOutput atts =
363  Application::GetUserInterface().GetOutputAttribute(parameter);
364 
365  FileName *file = new FileName(avgString);
366  QString path = file->path();
367  QString filename = file->baseName();
368  QString countString = path + "/" + filename + "-count";
369  SetStatCubes(avgString, countString, atts, nsamps, nlines, nbands);
370 
371  }
372 
373 
380  void ProcessPolygons::SetIntersectAlgorithm(const bool useCenter) {
381  m_useCenter = useCenter;
382  }
383 
384 
385 } /* end namespace isis*/
386 
Isis::BoxcarCachingAlgorithm
This algorithm is designed for applications that jump around between a couple of spots in the cube wi...
Definition: BoxcarCachingAlgorithm.h:32
Isis::FileName
File name manipulation and expansion.
Definition: FileName.h:100
Isis::CubeAttributeOutput
Manipulate and parse attributes of output cube filenames.
Definition: CubeAttribute.h:473
Isis::Brick
Buffer for containing a three dimensional section of an image.
Definition: Brick.h:45
Isis::FileName::baseName
QString baseName() const
Returns the name of the file without the path and without extensions.
Definition: FileName.cpp:145
Isis::Cube
IO Handler for Isis Cubes.
Definition: Cube.h:167
Isis::IException
Isis exception class.
Definition: IException.h:91
Isis::Null
const double Null
Value for an Isis Null pixel.
Definition: SpecialPixel.h:95
Isis::FileName::extension
QString extension() const
Returns the last extension of the file name.
Definition: FileName.cpp:178
std
Namespace for the standard library.
Isis::IsValidPixel
bool IsValidPixel(const double d)
Returns if the input pixel is valid.
Definition: SpecialPixel.h:223
Isis::Cube::open
void open(const QString &cfile, QString access="r")
This method will open an isis cube for reading or reading/writing.
Definition: Cube.cpp:627
Isis::FileName::path
QString path() const
Returns the path of the file name.
Definition: FileName.cpp:103
Isis
This is free and unencumbered software released into the public domain.
Definition: Apollo.h:16