10 #include "geos/operation/distance/DistanceOp.h"
11 #include "geos/util/IllegalArgumentException.h"
12 #include "geos/geom/Point.h"
13 #include "geos/opOverlay.h"
21 #include "QMessageBox"
36 ImageOverlapSet::ImageOverlapSet(
bool continueOnError) {
38 p_continueAfterError = continueOnError;
40 p_calculatedSoFar = -1;
41 p_threadedCalculate =
false;
52 ImageOverlapSet::~ImageOverlapSet() {
54 for (
int i = 0; i < Size(); i++) {
55 if (p_lonLatOverlaps[i])
delete p_lonLatOverlaps[i];
74 for (
int i = 0; i < sns.
size(); i++) {
81 QString msg =
"Unable to open cube for serial number [";
84 HandleError(error, &sns, msg);
93 geos::geom::MultiPolygon *tmp = PolygonTools::MakeMultiPolygon(poly->
Polys());
98 geos::geom::MultiPolygon *mp = NULL;
101 if(!tmp->isValid()) {
105 "] has an invalid footprint";
110 mp = PolygonTools::Despike(tmp);
113 if (tmp->isValid()) {
120 HandleError(e, &sns);
125 p_lonLatOverlapsMutex.lock();
126 p_lonLatOverlaps.push_back(CreateNewOverlap(sns.
serialNumber(i), mp));
127 p_lonLatOverlapsMutex.unlock();
142 DespikeLonLatOverlaps();
144 if (p_threadedCalculate) {
150 FindAllOverlaps(&sns);
170 void ImageOverlapSet::FindImageOverlaps(
SerialNumberList &boundaries, QString outputFile) {
173 if (!p_lonLatOverlaps.empty()) {
174 string msg =
"FindImageOverlaps(SerialNumberList&,QString) may not be called on " \
175 "an ImageOverlapSet which already contains overlaps.";
180 p_calculatedSoFar = -1;
182 p_snlist = &boundaries;
185 p_threadedCalculate =
true;
187 FindImageOverlaps(boundaries);
191 while (p_calculatedSoFar != p_lonLatOverlaps.size()) {
192 WriteImageOverlaps(outputFile);
196 if (p_calculatedSoFar != p_writtenSoFar)
197 WriteImageOverlaps(outputFile);
204 p_lonLatOverlaps.clear();
206 p_calculatedSoFar = -1;
207 p_threadedCalculate =
false;
300 void ImageOverlapSet::FindImageOverlaps(std::vector<QString> sns,
301 std::vector<geos::geom::MultiPolygon *> polygons) {
303 if (sns.size() != polygons.size()) {
304 string message =
"Invalid argument sizes. Sizes must match.";
309 for (
unsigned int i = 0; i < sns.size(); ++i) {
310 p_lonLatOverlapsMutex.lock();
311 p_lonLatOverlaps.push_back(CreateNewOverlap(sns[i], polygons[i]));
312 p_lonLatOverlapsMutex.unlock();
316 DespikeLonLatOverlaps();
328 void ImageOverlapSet::ReadImageOverlaps(
const QString &filename) {
330 QString file =
FileName(filename).expanded();
334 std::ifstream inStream;
335 QByteArray fileArray = file.toLatin1();
336 inStream.open(fileArray.constData(), fstream::in | fstream::binary);
338 while (!inStream.eof()) {
339 p_lonLatOverlapsMutex.lock();
341 p_lonLatOverlapsMutex.unlock();
347 p_lonLatOverlapsMutex.unlock();
348 QString msg =
"The overlap file [" + filename +
"] does not contain a "
349 "valid list of image overlaps";
353 p_lonLatOverlapsMutex.unlock();
354 QString msg =
"The overlap file [" + filename +
"] does not contain a "
355 "valid list of image overlaps";
378 bool ImageOverlapSet::SetPolygon(geos::geom::Geometry *poly,
383 bool success =
false;
384 geos::geom::MultiPolygon *multiPolygon = PolygonTools::MakeMultiPolygon(poly);
388 if (!multiPolygon->isValid() ||
389 (multiPolygon->getArea() < 1.0e-10 && !multiPolygon->isEmpty())) {
391 multiPolygon = Isis::globalFactory.createMultiPolygon();
394 if (position > p_lonLatOverlaps.size()) {
395 position = p_lonLatOverlaps.size();
399 if (!multiPolygon->isEmpty()) {
400 geos::geom::MultiPolygon *despiked = PolygonTools::Despike(multiPolygon);
402 multiPolygon = despiked;
408 if (multiPolygon->isValid() &&
409 (multiPolygon->isEmpty() || multiPolygon->getArea() > 1.0e-14)) {
411 p_lonLatOverlaps.at(position)->SetPolygon(multiPolygon);
414 AddSerialNumbers(p_lonLatOverlaps.at(position), sncopy);
417 else if (!multiPolygon->isEmpty()) {
424 AddSerialNumbers(imageOverlap, sncopy);
429 p_lonLatOverlapsMutex.lock();
430 p_lonLatOverlaps.insert(p_lonLatOverlaps.begin() + position, imageOverlap);
431 p_lonLatOverlapsMutex.unlock();
446 void ImageOverlapSet::WriteImageOverlaps(
const QString &filename) {
448 QString file =
FileName(filename).expanded();
450 bool noOverlaps =
false;
452 if (p_threadedCalculate) {
453 p_calculatePolygonMutex.lock();
458 std::ofstream outStream;
460 QByteArray fileArray = file.toLatin1();
461 if (p_writtenSoFar == 0) {
462 outStream.open(fileArray.constData(), fstream::out | fstream::trunc | fstream::binary);
465 outStream.open(fileArray.constData(), fstream::out | fstream::app | fstream::binary);
468 failed |= outStream.fail();
470 static bool overlapWritten =
false;
471 for (
int overlap = p_writtenSoFar; !failed && overlap <= p_calculatedSoFar; overlap++) {
473 p_lonLatOverlapsMutex.lock();
475 if (p_lonLatOverlaps.size() == 0) {
479 if (overlap < p_lonLatOverlaps.size() && p_lonLatOverlaps[overlap]) {
481 if (!p_lonLatOverlaps[overlap]->Polygon()->isEmpty()) {
483 if (overlapWritten) {
484 outStream << std::endl;
487 p_lonLatOverlaps[overlap]->Write(outStream);
488 overlapWritten =
true;
491 delete p_lonLatOverlaps[overlap];
492 p_lonLatOverlaps[overlap] = NULL;
496 p_lonLatOverlapsMutex.unlock();
499 failed |= outStream.fail();
502 failed |= outStream.fail();
512 if (p_calculatedSoFar == p_lonLatOverlaps.size()) {
513 if (p_threadedCalculate && !noOverlaps) {
514 p_calculatePolygonMutex.tryLock();
515 p_calculatePolygonMutex.unlock();
520 p_calculatePolygonMutex.tryLock();
521 p_calculatePolygonMutex.unlock();
522 QString msg =
"Unable to write the image overlap list to [" + filename +
"]";
525 else if (noOverlaps) {
526 p_calculatePolygonMutex.tryLock();
527 p_calculatePolygonMutex.unlock();
528 QString msg =
"No overlaps were found.";
542 bool foundOverlap =
false;
543 if (p_lonLatOverlaps.size() <= 1)
return;
546 p.
SetText(
"Calculating Image Overlaps");
550 geos::geom::MultiPolygon *emptyPolygon = Isis::globalFactory.createMultiPolygon();
553 for (
int outside = 0; outside < p_lonLatOverlaps.size() - 1; ++outside) {
554 p_calculatedSoFar = outside - 1;
557 if (p_calculatedSoFar % 10 == 0 && (!snlist || (p_lonLatOverlaps.size() > snlist->
size()))) {
558 if (p_threadedCalculate) {
559 p_calculatePolygonMutex.unlock();
565 for (
int inside = outside + 1; inside < p_lonLatOverlaps.size(); ++inside) {
567 if (p_lonLatOverlaps.at(outside)->HasAnySameSerialNumber(*p_lonLatOverlaps.at(inside)))
571 const geos::geom::MultiPolygon *poly1 = p_lonLatOverlaps.at(outside)->Polygon();
572 const geos::geom::MultiPolygon *poly2 = p_lonLatOverlaps.at(inside)->Polygon();
576 if (PolygonTools::Equal(poly1, poly2)) {
577 p_lonLatOverlapsMutex.lock();
578 AddSerialNumbers(p_lonLatOverlaps[outside], p_lonLatOverlaps[inside]);
579 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
580 p_lonLatOverlapsMutex.unlock();
586 if (poly2->isEmpty() || poly2->getArea() < 1.0e-14) {
587 p_lonLatOverlapsMutex.lock();
588 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
589 p_lonLatOverlapsMutex.unlock();
594 geos::geom::Geometry *intersected = NULL;
596 intersected = PolygonTools::Intersect(poly1, poly2);
600 QString error =
"Intersection of overlaps failed.";
605 double outsideArea = poly1->getArea();
606 double insideArea = poly2->getArea();
607 double areaRatio = std::min(outsideArea, insideArea) /
608 std::max(outsideArea, insideArea);
613 if (areaRatio < 0.1) {
614 if (poly1->getArea() > poly2->getArea()) {
615 error +=
" The first polygon will be removed.";
616 HandleError(e, snlist, error, inside, outside);
617 p_lonLatOverlapsMutex.lock();
618 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
619 p_lonLatOverlapsMutex.unlock();
623 error +=
" The second polygon will be removed.";
624 HandleError(e, snlist, error, inside, outside);
625 p_lonLatOverlapsMutex.lock();
626 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + outside);
627 p_lonLatOverlapsMutex.unlock();
632 error +=
" Both polygons will be removed to prevent the "
633 "possibility of double counted areas.";
634 HandleError(e, snlist, error, inside, outside);
635 p_lonLatOverlapsMutex.lock();
636 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
637 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + outside);
638 p_lonLatOverlapsMutex.unlock();
645 if (intersected->isEmpty() || intersected->getArea() < 1.0e-14) {
654 geos::geom::MultiPolygon *overlap = NULL;
656 overlap = PolygonTools::Despike(intersected);
662 if (!intersected->isValid()) {
666 HandleError(e, snlist,
"", inside, outside);
670 overlap = PolygonTools::MakeMultiPolygon(intersected);
676 catch (geos::util::GEOSException *exc) {
679 HandleError(exc, snlist,
"", inside, outside);
683 if (!overlap->isValid()) {
686 HandleError(snlist,
"Intersection produced invalid overlap area", inside, outside);
691 if (overlap->isEmpty() || overlap->getArea() < 1.0e-14) {
698 if (PolygonTools::Equal(poly1, overlap)) {
699 geos::geom::Geometry *tmpGeom = NULL;
701 tmpGeom = PolygonTools::Difference(poly2, poly1);
704 HandleError(e, snlist,
"Differencing overlap polygons failed."
705 "The first polygon will be removed.", inside, outside);
709 p_lonLatOverlapsMutex.lock();
710 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + outside);
711 p_lonLatOverlapsMutex.unlock();
715 if (SetPolygon(tmpGeom, inside) &&
716 SetPolygon(overlap, outside, p_lonLatOverlaps[inside]))
720 else if (PolygonTools::Equal(poly2, overlap)) {
721 geos::geom::Geometry *tmpGeom = NULL;
723 tmpGeom = PolygonTools::Difference(poly1, poly2);
726 HandleError(e, snlist,
"Differencing overlap polygons failed."
727 "The second polygon will be removed.", inside, outside);
730 p_lonLatOverlapsMutex.lock();
731 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
732 p_lonLatOverlapsMutex.unlock();
736 if (SetPolygon(tmpGeom, outside) &&
737 SetPolygon(overlap, inside, p_lonLatOverlaps[outside]))
743 geos::geom::Geometry *tmpGeom = NULL;
745 tmpGeom = PolygonTools::Difference(poly1, overlap);
754 if (tmpGeom == NULL) {
755 tmpGeom = PolygonTools::Difference(poly1, poly2);
760 HandleError(e, snlist,
"Differencing overlap polygons failed", inside, outside);
764 if (!SetPolygon(tmpGeom, outside)) {
765 if (SetPolygon(Isis::globalFactory.createMultiPolygon(), outside))
769 int oldSize = p_lonLatOverlaps.size();
770 if (SetPolygon(overlap, inside + 1, p_lonLatOverlaps[outside],
true)) {
771 int newSize = p_lonLatOverlaps.size();
772 int newSteps = newSize - oldSize;
775 if (newSize != oldSize) inside++;
781 HandleError(e, snlist,
"Unable to find overlap.", inside, outside);
784 catch (geos::util::IllegalArgumentException *ill) {
785 HandleError(NULL, snlist,
"Unable to find overlap", inside, outside);
787 catch (geos::util::GEOSException *exc) {
788 HandleError(exc, snlist,
"Unable to find overlap", inside, outside);
791 HandleError(snlist,
"Unknown Error: Unable to find overlap", inside, outside);
798 p_calculatedSoFar = p_lonLatOverlaps.size();
802 if (foundOverlap ==
false) {
803 p_lonLatOverlapsMutex.lock();
804 p_lonLatOverlaps.clear();
805 p_lonLatOverlapsMutex.unlock();
811 p_calculatePolygonMutex.tryLock();
812 p_calculatePolygonMutex.unlock();
826 for (
int i = 0; i < from->Size(); i++) {
827 QString s = (*from)[i];
840 geos::geom::MultiPolygon *latLonPolygon) {
857 std::vector<ImageOverlap *> ImageOverlapSet::operator[](QString serialNumber) {
859 vector<ImageOverlap *> matches;
863 for (
int ov = 0; ov < p_lonLatOverlaps.size(); ++ov) {
864 for (
int sn = 0; sn < p_lonLatOverlaps[ov]->Size(); ++sn) {
865 if ((*p_lonLatOverlaps[ov])[sn] == serialNumber) {
866 matches.push_back(p_lonLatOverlaps[ov]);
888 int overlap1,
int overlap2) {
892 if (overlap1 >= 0 && overlap1 < p_lonLatOverlaps.size()) {
893 PvlKeyword serialNumbers(
"PolySerialNumbers");
897 for (
int i = 0; i < p_lonLatOverlaps.at(overlap1)->Size(); i++) {
898 serialNumbers += (*p_lonLatOverlaps.at(overlap1))[i];
900 if (snlist != NULL) {
901 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap1))[i]);
904 polygon += p_lonLatOverlaps.at(overlap1)->Polygon()->toString().c_str();
906 err += serialNumbers;
908 if (filename.
size() != 0) {
915 if (overlap2 >= 0 && overlap1 < p_lonLatOverlaps.size() &&
916 overlap2 < p_lonLatOverlaps.size()) {
917 PvlKeyword serialNumbers(
"PolySerialNumbers");
921 for (
int i = 0; i < p_lonLatOverlaps.at(overlap2)->Size(); i++) {
922 serialNumbers += (*p_lonLatOverlaps.at(overlap2))[i];
924 if (snlist != NULL) {
925 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap2))[i]);
928 polygon += p_lonLatOverlaps.at(overlap2)->Polygon()->toString().c_str();
930 err += serialNumbers;
932 if (filename.
size() != 0) {
941 if (!msg.isEmpty()) {
945 p_errorLog.push_back(err);
947 if (!p_continueAfterError)
throw;
961 void ImageOverlapSet::HandleError(geos::util::GEOSException *exc,
964 int overlap1,
int overlap2) {
968 if (overlap1 >= 0 && overlap1 < p_lonLatOverlaps.size()) {
969 PvlKeyword serialNumbers(
"PolySerialNumbers");
972 for (
int i = 0; i < p_lonLatOverlaps.at(overlap1)->Size(); i++) {
973 serialNumbers += (*p_lonLatOverlaps.at(overlap1))[i];
975 if (snlist != NULL) {
976 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap1))[i]);
980 err += serialNumbers;
982 if (filename.
size() != 0) {
987 if (overlap2 >= 0 && overlap1 < p_lonLatOverlaps.size() &&
988 overlap2 < p_lonLatOverlaps.size()) {
989 PvlKeyword serialNumbers(
"PolySerialNumbers");
992 for (
int i = 0; i < p_lonLatOverlaps.at(overlap2)->Size(); i++) {
993 serialNumbers += (*p_lonLatOverlaps.at(overlap2))[i];
995 if (snlist != NULL) {
996 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap2))[i]);
1000 err += serialNumbers;
1002 if (filename.
size() != 0) {
1009 if (!msg.isEmpty()) {
1013 p_errorLog.push_back(err);
1017 if (!p_continueAfterError) {
1034 int overlap1,
int overlap2) {
1038 if (overlap1 >= 0 && overlap1 < p_lonLatOverlaps.size()) {
1039 PvlKeyword serialNumbers(
"PolySerialNumbers");
1042 for (
int i = 0; i < p_lonLatOverlaps.at(overlap1)->Size(); i++) {
1043 serialNumbers += (*p_lonLatOverlaps.at(overlap1))[i];
1045 if (snlist != NULL) {
1046 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap1))[i]);
1050 err += serialNumbers;
1052 if (filename.
size() != 0) {
1057 if (overlap2 >= 0 && overlap1 < p_lonLatOverlaps.size() &&
1058 overlap2 < p_lonLatOverlaps.size()) {
1059 PvlKeyword serialNumbers(
"PolySerialNumbers");
1062 for (
int i = 0; i < p_lonLatOverlaps.at(overlap2)->Size(); i++) {
1063 serialNumbers += (*p_lonLatOverlaps.at(overlap2))[i];
1065 if (snlist != NULL) {
1066 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap2))[i]);
1070 err += serialNumbers;
1072 if (filename.
size() != 0) {
1079 p_errorLog.push_back(err);
1081 if (!p_continueAfterError) {
1092 void ImageOverlapSet::DespikeLonLatOverlaps() {
1094 for (
int i = 0; i < Size(); i++) {
1096 p_lonLatOverlaps[i]->SetPolygon(PolygonTools::Despike(p_lonLatOverlaps[i]->Polygon()));
int size() const
Returns the number of values stored in this keyword.
void SetMaximumSteps(const int steps)
This sets the maximum number of steps in the process.
File name manipulation and expansion.
int size() const
How many serial number / filename combos are in the list.
QString serialNumber(const QString &filename)
Return a serial number given a filename.
void Add(QString &sn)
This method will add a new serial number to the list of serial numbers alread associated with the ove...
virtual void SetPolygon(const geos::geom::MultiPolygon &polygon)
This method will replace the existing polygon that defines the overlap with a new one...
void read(Blob &blob) const
This method will read data from the specified Blob object.
Create cube polygons, read/write polygons to blobs.
void CheckStatus()
Checks and updates the status.
const char * what() const
Returns a string representation of this exception in its current state.
QString fileName(const QString &sn)
Return a filename given a serial number.
Program progress reporter.
void AddSteps(const int steps)
If the initial step size was a guess, it can be modified using this method.
void SetText(const QString &text)
Changes the value of the text string reported just before 0% processed.
Contains multiple PvlContainers.
#define _FILEINFO_
Macro for the filename and line number.
A single keyword-value pair.
void close(bool remove=false)
Closes the cube and updates the labels.
void open(const QString &cfile, QString access="r")
This method will open an isis cube for reading or reading/writing.
Individual overlap container.
geos::geom::MultiPolygon * Polys()
Return a geos Multipolygon.
Serial Number list generator.
IO Handler for Isis Cubes.