16 #include "geos/operation/distance/DistanceOp.h"
17 #include "geos/util/IllegalArgumentException.h"
18 #include "geos/geom/Point.h"
19 #include "geos/opOverlay.h"
20 #include "IException.h"
21 #include "ImageOverlapSet.h"
22 #include "ImagePolygon.h"
23 #include "PolygonTools.h"
25 #include "SerialNumberList.h"
27 #include "QMessageBox"
42 ImageOverlapSet::ImageOverlapSet(
bool continueOnError,
bool useThread) {
44 p_continueAfterError = continueOnError;
46 p_calculatedSoFar = -1;
47 p_threadedCalculate = useThread;
58 ImageOverlapSet::~ImageOverlapSet() {
60 for (
int i = 0; i < Size(); i++) {
61 if (p_lonLatOverlaps[i])
delete p_lonLatOverlaps[i];
80 for (
int i = 0; i < sns.
size(); i++) {
87 QString msg =
"Unable to open cube for serial number [";
90 HandleError(error, &sns, msg);
98 geos::geom::MultiPolygon *tmp = PolygonTools::MakeMultiPolygon(poly.
Polys());
100 geos::geom::MultiPolygon *mp = NULL;
103 if(!tmp->isValid()) {
107 "] has an invalid footprint";
108 throw IException(IException::Programmer, msg, _FILEINFO_);
112 mp = PolygonTools::Despike(tmp);
115 if (tmp->isValid()) {
122 HandleError(e, &sns);
127 p_lonLatOverlapsMutex.lock();
128 p_lonLatOverlaps.push_back(CreateNewOverlap(sns.
serialNumber(i), mp));
129 p_lonLatOverlapsMutex.unlock();
144 DespikeLonLatOverlaps();
146 if (p_threadedCalculate) {
152 FindAllOverlaps(&sns);
172 void ImageOverlapSet::FindImageOverlaps(
SerialNumberList &boundaries, QString outputFile) {
175 if (!p_lonLatOverlaps.empty()) {
176 string msg =
"FindImageOverlaps(SerialNumberList&,QString) may not be called on " \
177 "an ImageOverlapSet which already contains overlaps.";
178 throw IException(IException::Programmer, msg, _FILEINFO_);
182 p_calculatedSoFar = -1;
184 p_snlist = &boundaries;
186 FindImageOverlaps(boundaries);
190 while (p_calculatedSoFar != p_lonLatOverlaps.size()) {
191 WriteImageOverlaps(outputFile);
195 if (p_calculatedSoFar != p_writtenSoFar)
196 WriteImageOverlaps(outputFile);
203 p_lonLatOverlaps.clear();
205 p_calculatedSoFar = -1;
206 p_threadedCalculate =
false;
299 void ImageOverlapSet::FindImageOverlaps(std::vector<QString> sns,
300 std::vector<geos::geom::MultiPolygon *> polygons) {
302 if (sns.size() != polygons.size()) {
303 string message =
"Invalid argument sizes. Sizes must match.";
304 throw IException(IException::Programmer, message, _FILEINFO_);
308 for (
unsigned int i = 0; i < sns.size(); ++i) {
309 p_lonLatOverlapsMutex.lock();
310 p_lonLatOverlaps.push_back(CreateNewOverlap(sns[i], polygons[i]));
311 p_lonLatOverlapsMutex.unlock();
315 DespikeLonLatOverlaps();
327 void ImageOverlapSet::ReadImageOverlaps(
const QString &filename) {
333 std::ifstream inStream;
334 QByteArray fileArray = file.toLatin1();
335 inStream.open(fileArray.constData(), fstream::in | fstream::binary);
337 while (!inStream.eof()) {
338 p_lonLatOverlapsMutex.lock();
340 p_lonLatOverlapsMutex.unlock();
346 p_lonLatOverlapsMutex.unlock();
347 QString msg =
"The overlap file [" + filename +
"] does not contain a "
348 "valid list of image overlaps";
349 throw IException(e, IException::Unknown, msg, _FILEINFO_);
352 p_lonLatOverlapsMutex.unlock();
353 QString msg =
"The overlap file [" + filename +
"] does not contain a "
354 "valid list of image overlaps";
355 throw IException(IException::Unknown, msg, _FILEINFO_);
377 bool ImageOverlapSet::SetPolygon(geos::geom::Geometry *poly,
382 bool success =
false;
383 geos::geom::MultiPolygon *multiPolygon = PolygonTools::MakeMultiPolygon(poly);
387 if (!multiPolygon->isValid() ||
388 (multiPolygon->getArea() < 1.0e-10 && !multiPolygon->isEmpty())) {
390 multiPolygon = Isis::globalFactory->createMultiPolygon();
393 if (position > p_lonLatOverlaps.size()) {
394 position = p_lonLatOverlaps.size();
398 if (!multiPolygon->isEmpty()) {
399 geos::geom::MultiPolygon *despiked = PolygonTools::Despike(multiPolygon);
401 multiPolygon = despiked;
407 if (multiPolygon->isValid() &&
408 (multiPolygon->isEmpty() || multiPolygon->getArea() > 1.0e-14)) {
410 p_lonLatOverlaps.at(position)->SetPolygon(multiPolygon);
413 AddSerialNumbers(p_lonLatOverlaps.at(position), sncopy);
416 else if (!multiPolygon->isEmpty()) {
423 AddSerialNumbers(imageOverlap, sncopy);
428 p_lonLatOverlapsMutex.lock();
429 p_lonLatOverlaps.insert(p_lonLatOverlaps.begin() + position, imageOverlap);
430 p_lonLatOverlapsMutex.unlock();
445 void ImageOverlapSet::WriteImageOverlaps(
const QString &filename) {
449 bool noOverlaps =
false;
451 if (p_threadedCalculate) {
452 p_calculatePolygonMutex.lock();
456 std::ofstream outStream;
459 QByteArray fileArray = file.toLatin1();
460 if (p_writtenSoFar == 0) {
461 outStream.open(fileArray.constData(), fstream::out | fstream::trunc | fstream::binary);
464 outStream.open(fileArray.constData(), fstream::out | fstream::app | fstream::binary);
467 failed |= outStream.fail();
469 for (
int overlap = p_writtenSoFar; !failed && overlap <= p_calculatedSoFar; overlap++) {
471 p_lonLatOverlapsMutex.lock();
473 if (p_lonLatOverlaps.size() == 0) {
477 if (overlap < p_lonLatOverlaps.size() && p_lonLatOverlaps[overlap]) {
479 if (!p_lonLatOverlaps[overlap]->Polygon()->isEmpty()) {
481 if (p_writtenSoFar) {
482 outStream << std::endl;
485 p_lonLatOverlaps[overlap]->Write(outStream);
488 delete p_lonLatOverlaps[overlap];
489 p_lonLatOverlaps[overlap] = NULL;
493 p_lonLatOverlapsMutex.unlock();
496 failed |= outStream.fail();
499 failed |= outStream.fail();
510 if (p_calculatedSoFar == p_lonLatOverlaps.size()) {
511 if (p_threadedCalculate && !noOverlaps) {
512 p_calculatePolygonMutex.tryLock();
513 p_calculatePolygonMutex.unlock();
518 p_calculatePolygonMutex.tryLock();
519 p_calculatePolygonMutex.unlock();
520 QString msg =
"Unable to write the image overlap list to [" + filename +
"]";
521 throw IException(IException::Io, msg, _FILEINFO_);
523 else if (noOverlaps) {
524 p_calculatePolygonMutex.tryLock();
525 p_calculatePolygonMutex.unlock();
526 QString msg =
"No overlaps were found.";
527 throw IException(IException::User, msg, _FILEINFO_);
540 bool foundOverlap =
false;
541 if (p_lonLatOverlaps.size() <= 1)
return;
544 p.
SetText(
"Calculating Image Overlaps");
548 geos::geom::MultiPolygon *emptyPolygon = Isis::globalFactory->createMultiPolygon();
551 for (
int outside = 0; outside < p_lonLatOverlaps.size() - 1; ++outside) {
552 p_calculatedSoFar = outside - 1;
555 if (p_calculatedSoFar % 10 == 0 && (!snlist || (p_lonLatOverlaps.size() > snlist->
size()))) {
556 if (p_threadedCalculate) {
557 p_calculatePolygonMutex.tryLock();
558 p_calculatePolygonMutex.unlock();
564 for (
int inside = outside + 1; inside < p_lonLatOverlaps.size(); ++inside) {
566 if (p_lonLatOverlaps.at(outside)->HasAnySameSerialNumber(*p_lonLatOverlaps.at(inside)))
570 const geos::geom::MultiPolygon *poly1 = p_lonLatOverlaps.at(outside)->Polygon();
571 const geos::geom::MultiPolygon *poly2 = p_lonLatOverlaps.at(inside)->Polygon();
575 if (PolygonTools::Equal(poly1, poly2)) {
576 p_lonLatOverlapsMutex.lock();
577 AddSerialNumbers(p_lonLatOverlaps[outside], p_lonLatOverlaps[inside]);
578 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
579 p_lonLatOverlapsMutex.unlock();
585 if (poly2->isEmpty() || poly2->getArea() < 1.0e-14) {
586 p_lonLatOverlapsMutex.lock();
587 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
588 p_lonLatOverlapsMutex.unlock();
593 geos::geom::Geometry *intersected = NULL;
595 intersected = PolygonTools::Intersect(poly1, poly2);
599 QString error =
"Intersection of overlaps failed.";
604 double outsideArea = poly1->getArea();
605 double insideArea = poly2->getArea();
606 double areaRatio = std::min(outsideArea, insideArea) /
607 std::max(outsideArea, insideArea);
612 if (areaRatio < 0.1) {
613 if (poly1->getArea() > poly2->getArea()) {
614 error +=
" The first polygon will be removed.";
615 HandleError(e, snlist, error, inside, outside);
616 p_lonLatOverlapsMutex.lock();
617 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
618 p_lonLatOverlapsMutex.unlock();
622 error +=
" The second polygon will be removed.";
623 HandleError(e, snlist, error, inside, outside);
624 p_lonLatOverlapsMutex.lock();
625 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + outside);
626 p_lonLatOverlapsMutex.unlock();
631 error +=
" Both polygons will be removed to prevent the "
632 "possibility of double counted areas.";
633 HandleError(e, snlist, error, inside, outside);
634 p_lonLatOverlapsMutex.lock();
635 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
636 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + outside);
637 p_lonLatOverlapsMutex.unlock();
644 if (intersected->isEmpty() || intersected->getArea() < 1.0e-14) {
653 geos::geom::MultiPolygon *overlap = NULL;
655 overlap = PolygonTools::Despike(intersected);
661 if (!intersected->isValid()) {
665 HandleError(e, snlist,
"", inside, outside);
669 overlap = PolygonTools::MakeMultiPolygon(intersected);
675 catch (geos::util::GEOSException *exc) {
678 HandleError(exc, snlist,
"", inside, outside);
682 if (!overlap->isValid()) {
685 HandleError(snlist,
"Intersection produced invalid overlap area", inside, outside);
690 if (overlap->isEmpty() || overlap->getArea() < 1.0e-14) {
697 if (PolygonTools::Equal(poly1, overlap)) {
698 geos::geom::Geometry *tmpGeom = NULL;
700 tmpGeom = PolygonTools::Difference(poly2, poly1);
703 HandleError(e, snlist,
"Differencing overlap polygons failed."
704 "The first polygon will be removed.", inside, outside);
708 p_lonLatOverlapsMutex.lock();
709 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + outside);
710 p_lonLatOverlapsMutex.unlock();
714 if (SetPolygon(tmpGeom, inside) &&
715 SetPolygon(overlap, outside, p_lonLatOverlaps[inside]))
719 else if (PolygonTools::Equal(poly2, overlap)) {
720 geos::geom::Geometry *tmpGeom = NULL;
722 tmpGeom = PolygonTools::Difference(poly1, poly2);
725 HandleError(e, snlist,
"Differencing overlap polygons failed."
726 "The second polygon will be removed.", inside, outside);
729 p_lonLatOverlapsMutex.lock();
730 p_lonLatOverlaps.erase(p_lonLatOverlaps.begin() + inside);
731 p_lonLatOverlapsMutex.unlock();
735 if (SetPolygon(tmpGeom, outside) &&
736 SetPolygon(overlap, inside, p_lonLatOverlaps[outside]))
742 geos::geom::Geometry *tmpGeom = NULL;
744 tmpGeom = PolygonTools::Difference(poly1, overlap);
753 if (tmpGeom == NULL) {
754 tmpGeom = PolygonTools::Difference(poly1, poly2);
759 HandleError(e, snlist,
"Differencing overlap polygons failed", inside, outside);
763 if (!SetPolygon(tmpGeom, outside)) {
764 if (SetPolygon(Isis::globalFactory->createMultiPolygon(), outside))
768 int oldSize = p_lonLatOverlaps.size();
769 if (SetPolygon(overlap, inside + 1, p_lonLatOverlaps[outside],
true)) {
770 int newSize = p_lonLatOverlaps.size();
771 int newSteps = newSize - oldSize;
774 if (newSize != oldSize) inside++;
780 HandleError(e, snlist,
"Unable to find overlap.", inside, outside);
783 catch (geos::util::IllegalArgumentException *ill) {
784 HandleError(NULL, snlist,
"Unable to find overlap", inside, outside);
786 catch (geos::util::GEOSException *exc) {
787 HandleError(exc, snlist,
"Unable to find overlap", inside, outside);
790 HandleError(snlist,
"Unknown Error: Unable to find overlap", inside, outside);
797 p_calculatedSoFar = p_lonLatOverlaps.size();
801 if (foundOverlap ==
false) {
802 p_lonLatOverlapsMutex.lock();
803 p_lonLatOverlaps.clear();
804 p_lonLatOverlapsMutex.unlock();
810 p_calculatePolygonMutex.tryLock();
811 p_calculatePolygonMutex.unlock();
825 for (
int i = 0; i < from->Size(); i++) {
826 QString s = (*from)[i];
839 geos::geom::MultiPolygon *latLonPolygon) {
856 std::vector<ImageOverlap *> ImageOverlapSet::operator[](QString serialNumber) {
858 vector<ImageOverlap *> matches;
862 for (
int ov = 0; ov < p_lonLatOverlaps.size(); ++ov) {
863 for (
int sn = 0; sn < p_lonLatOverlaps[ov]->Size(); ++sn) {
864 if ((*p_lonLatOverlaps[ov])[sn] == serialNumber) {
865 matches.push_back(p_lonLatOverlaps[ov]);
887 int overlap1,
int overlap2) {
891 if (overlap1 >= 0 && overlap1 < p_lonLatOverlaps.size()) {
892 PvlKeyword serialNumbers(
"PolySerialNumbers");
896 for (
int i = 0; i < p_lonLatOverlaps.at(overlap1)->Size(); i++) {
897 serialNumbers += (*p_lonLatOverlaps.at(overlap1))[i];
899 if (snlist != NULL) {
900 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap1))[i]);
903 polygon += p_lonLatOverlaps.at(overlap1)->Polygon()->toString().c_str();
905 err += serialNumbers;
907 if (filename.
size() != 0) {
914 if (overlap2 >= 0 && overlap1 < p_lonLatOverlaps.size() &&
915 overlap2 < p_lonLatOverlaps.size()) {
916 PvlKeyword serialNumbers(
"PolySerialNumbers");
920 for (
int i = 0; i < p_lonLatOverlaps.at(overlap2)->Size(); i++) {
921 serialNumbers += (*p_lonLatOverlaps.at(overlap2))[i];
923 if (snlist != NULL) {
924 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap2))[i]);
927 polygon += p_lonLatOverlaps.at(overlap2)->Polygon()->toString().c_str();
929 err += serialNumbers;
931 if (filename.
size() != 0) {
940 if (!msg.isEmpty()) {
944 p_errorLog.push_back(err);
946 if (!p_continueAfterError)
throw;
960 void ImageOverlapSet::HandleError(geos::util::GEOSException *exc,
963 int overlap1,
int overlap2) {
967 if (overlap1 >= 0 && overlap1 < p_lonLatOverlaps.size()) {
968 PvlKeyword serialNumbers(
"PolySerialNumbers");
971 for (
int i = 0; i < p_lonLatOverlaps.at(overlap1)->Size(); i++) {
972 serialNumbers += (*p_lonLatOverlaps.at(overlap1))[i];
974 if (snlist != NULL) {
975 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap1))[i]);
979 err += serialNumbers;
981 if (filename.
size() != 0) {
986 if (overlap2 >= 0 && overlap1 < p_lonLatOverlaps.size() &&
987 overlap2 < p_lonLatOverlaps.size()) {
988 PvlKeyword serialNumbers(
"PolySerialNumbers");
991 for (
int i = 0; i < p_lonLatOverlaps.at(overlap2)->Size(); i++) {
992 serialNumbers += (*p_lonLatOverlaps.at(overlap2))[i];
994 if (snlist != NULL) {
995 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap2))[i]);
999 err += serialNumbers;
1001 if (filename.
size() != 0) {
1008 if (!msg.isEmpty()) {
1012 p_errorLog.push_back(err);
1016 if (!p_continueAfterError) {
1017 throw IException(IException::Programmer, err[
"Description"][0], _FILEINFO_);
1033 int overlap1,
int overlap2) {
1037 if (overlap1 >= 0 && overlap1 < p_lonLatOverlaps.size()) {
1038 PvlKeyword serialNumbers(
"PolySerialNumbers");
1041 for (
int i = 0; i < p_lonLatOverlaps.at(overlap1)->Size(); i++) {
1042 serialNumbers += (*p_lonLatOverlaps.at(overlap1))[i];
1044 if (snlist != NULL) {
1045 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap1))[i]);
1049 err += serialNumbers;
1051 if (filename.
size() != 0) {
1056 if (overlap2 >= 0 && overlap1 < p_lonLatOverlaps.size() &&
1057 overlap2 < p_lonLatOverlaps.size()) {
1058 PvlKeyword serialNumbers(
"PolySerialNumbers");
1061 for (
int i = 0; i < p_lonLatOverlaps.at(overlap2)->Size(); i++) {
1062 serialNumbers += (*p_lonLatOverlaps.at(overlap2))[i];
1064 if (snlist != NULL) {
1065 filename += snlist->
fileName((*p_lonLatOverlaps.at(overlap2))[i]);
1069 err += serialNumbers;
1071 if (filename.
size() != 0) {
1078 p_errorLog.push_back(err);
1080 if (!p_continueAfterError) {
1081 throw IException(IException::Programmer, err[
"Description"][0], _FILEINFO_);
1091 void ImageOverlapSet::DespikeLonLatOverlaps() {
1093 for (
int i = 0; i < Size(); i++) {
1095 p_lonLatOverlaps[i]->SetPolygon(PolygonTools::Despike(p_lonLatOverlaps[i]->Polygon()));