diff --git a/src/engine/SpatialJoin.cpp b/src/engine/SpatialJoin.cpp index beb09ffe6..4bf1eb372 100644 --- a/src/engine/SpatialJoin.cpp +++ b/src/engine/SpatialJoin.cpp @@ -209,7 +209,7 @@ size_t SpatialJoin::getCostEstimate() { if (childLeft_ && childRight_) { size_t inputEstimate = childLeft_->getSizeEstimate() * childRight_->getSizeEstimate(); - if (useBaselineAlgorithm_) { + if (algorithm_ == Algorithm::Baseline) { return inputEstimate * inputEstimate; } else { // Let n be the size of the left table and m the size of the right table. @@ -337,10 +337,21 @@ PreparedSpatialJoinParams SpatialJoin::prepareJoin() const { Result SpatialJoin::computeResult([[maybe_unused]] bool requestLaziness) { SpatialJoinAlgorithms algorithms{_executionContext, prepareJoin(), addDistToResult_, config_}; - if (useBaselineAlgorithm_) { + if (algorithm_ == Algorithm::Baseline) { return algorithms.BaselineAlgorithm(); - } else { + } else if (algorithm_ == Algorithm::S2Geometry) { return algorithms.S2geometryAlgorithm(); + } else if (algorithm_ == Algorithm::BoundingBox) { + // as the BoundingBoxAlgorithms only works for max distance and not for + // nearest neighbors, S2geometry gets called as a backup, if the query is + // asking for the nearest neighbors + if (std::get_if(&config_)) { + return algorithms.BoundingBoxAlgorithm(); + } else { + return algorithms.S2geometryAlgorithm(); + } + } else { + AD_THROW("choose a valid Algorithm"); } } diff --git a/src/engine/SpatialJoin.h b/src/engine/SpatialJoin.h index b62978cf9..8c6768223 100644 --- a/src/engine/SpatialJoin.h +++ b/src/engine/SpatialJoin.h @@ -96,14 +96,20 @@ class SpatialJoin : public Operation { // purposes std::optional getMaxResults() const; - void selectAlgorithm(bool useBaselineAlgorithm) { - useBaselineAlgorithm_ = useBaselineAlgorithm; - } + // options which can be used for the algorithm, which calculates the result + enum Algorithm { Baseline, S2Geometry, BoundingBox }; + + void selectAlgorithm(Algorithm algorithm) { algorithm_ = algorithm; } std::pair onlyForTestingGetConfig() const { return std::pair{getMaxDist().value_or(-1), getMaxResults().value_or(-1)}; } + std::variant + onlyForTestingGetRealConfig() const { + return config_; + } + std::shared_ptr onlyForTestingGetLeftChild() const { return childLeft_; } @@ -135,5 +141,5 @@ class SpatialJoin : public Operation { // between the two objects bool addDistToResult_ = true; const string nameDistanceInternal_ = "?distOfTheTwoObjectsAddedInternally"; - bool useBaselineAlgorithm_ = false; + Algorithm algorithm_ = BoundingBox; }; diff --git a/src/engine/SpatialJoinAlgorithms.cpp b/src/engine/SpatialJoinAlgorithms.cpp index dd7a9b1c7..2935d9000 100644 --- a/src/engine/SpatialJoinAlgorithms.cpp +++ b/src/engine/SpatialJoinAlgorithms.cpp @@ -223,3 +223,267 @@ Result SpatialJoinAlgorithms::S2geometryAlgorithm() { return Result(std::move(result), std::vector{}, Result::getMergedLocalVocab(*resultLeft, *resultRight)); } + +// ____________________________________________________________________________ +std::vector SpatialJoinAlgorithms::computeBoundingBox( + const point& startPoint) { + const auto [idTableLeft, resultLeft, idTableRight, resultRight, leftJoinCol, + rightJoinCol, numColumns, maxDist, maxResults] = params_; + if (!maxDist.has_value()) { + AD_THROW("Max distance must have a value for this operation"); + } + // haversine function + auto haversine = [](double theta) { return (1 - std::cos(theta)) / 2; }; + + // inverse haversine function + auto archaversine = [](double theta) { return std::acos(1 - 2 * theta); }; + + // safety buffer for numerical inaccuracies + double maxDistInMetersBuffer = static_cast(maxDist.value()); + if (maxDist.value() < 10) { + maxDistInMetersBuffer = 10; + } else if (maxDist.value() < std::numeric_limits::max() / 1.02) { + maxDistInMetersBuffer = 1.01 * maxDist.value(); + } else { + maxDistInMetersBuffer = std::numeric_limits::max(); + } + + // for large distances, where the lower calculation would just result in + // a single bounding box for the whole planet, do an optimized version + if (maxDist.value() > circumferenceMax_ / 4.0 && + maxDist.value() < circumferenceMax_ / 2.01) { + return computeAntiBoundingBox(startPoint); + } + + // compute latitude bound + double upperLatBound = + startPoint.get<1>() + maxDistInMetersBuffer * (360 / circumferenceMax_); + double lowerLatBound = + startPoint.get<1>() - maxDistInMetersBuffer * (360 / circumferenceMax_); + bool poleReached = false; + // test for "overflows" + if (lowerLatBound <= -90) { + lowerLatBound = -90; + poleReached = true; // south pole reached + } + if (upperLatBound >= 90) { + upperLatBound = 90; + poleReached = true; // north pole reached + } + if (poleReached) { + return {box(point(-180.0f, lowerLatBound), point(180.0f, upperLatBound))}; + } + + // compute longitude bound. For an explanation of the calculation and the + // naming convention see my master thesis + double alpha = maxDistInMetersBuffer / radius_; + double gamma = + (90 - std::abs(startPoint.get<1>())) * (2 * std::numbers::pi / 360); + double beta = std::acos((std::cos(gamma) / std::cos(alpha))); + double delta = 0; + if (maxDistInMetersBuffer > circumferenceMax_ / 20) { + // use law of cosines + delta = std::acos((std::cos(alpha) - std::cos(gamma) * std::cos(beta)) / + (std::sin(gamma) * std::sin(beta))); + } else { + // use law of haversines for numerical stability + delta = archaversine((haversine(alpha - haversine(gamma - beta))) / + (std::sin(gamma) * std::sin(beta))); + } + double lonRange = delta * 360 / (2 * std::numbers::pi); + double leftLonBound = startPoint.get<0>() - lonRange; + double rightLonBound = startPoint.get<0>() + lonRange; + // test for "overflows" and create two bounding boxes if necessary + if (leftLonBound < -180) { + box box1 = + box(point(-180, lowerLatBound), point(rightLonBound, upperLatBound)); + box box2 = box(point(leftLonBound + 360, lowerLatBound), + point(180, upperLatBound)); + return {box1, box2}; + } else if (rightLonBound > 180) { + box box1 = + box(point(leftLonBound, lowerLatBound), point(180, upperLatBound)); + box box2 = box(point(-180, lowerLatBound), + point(rightLonBound - 360, upperLatBound)); + return {box1, box2}; + } + // default case, when no bound has an "overflow" + return {box(point(leftLonBound, lowerLatBound), + point(rightLonBound, upperLatBound))}; +} + +// ____________________________________________________________________________ +std::vector SpatialJoinAlgorithms::computeAntiBoundingBox( + const point& startPoint) { + const auto [idTableLeft, resultLeft, idTableRight, resultRight, leftJoinCol, + rightJoinCol, numColumns, maxDist, maxResults] = params_; + if (!maxDist.has_value()) { + AD_THROW("Max distance must have a value for this operation"); + } + // point on the opposite side of the globe + point antiPoint(startPoint.get<0>() + 180, startPoint.get<1>() * -1); + if (antiPoint.get<0>() > 180) { + antiPoint.set<0>(antiPoint.get<0>() - 360); + } + // for an explanation of the formula see the master thesis. Divide by two two + // only consider the distance from the point to the antiPoint, subtract + // maxDist and a safety margine from that + double antiDist = + (circumferenceMin_ / 2.0) - maxDist.value() * 1.01; // safety margin + // use the bigger circumference as an additional safety margin, use 2.01 + // instead of 2.0 because of rounding inaccuracies in floating point + // operations + double distToAntiPoint = (360 / circumferenceMax_) * (antiDist / 2.01); + double upperBound = antiPoint.get<1>() + distToAntiPoint; + double lowerBound = antiPoint.get<1>() - distToAntiPoint; + double leftBound = antiPoint.get<0>() - distToAntiPoint; + double rightBound = antiPoint.get<0>() + distToAntiPoint; + bool northPoleTouched = false; + bool southPoleTouched = false; + bool boxCrosses180Longitude = false; // if the 180 to -180 line is touched + // if a pole is crossed, ignore the part after the crossing + if (upperBound > 90) { + upperBound = 90; + northPoleTouched = true; + } + if (lowerBound < -90) { + lowerBound = -90; + southPoleTouched = true; + } + if (leftBound < -180) { + leftBound += 360; + } + if (rightBound > 180) { + rightBound -= 360; + } + if (rightBound < leftBound) { + boxCrosses180Longitude = true; + } + // compute bounding boxes using the anti bounding box from above + std::vector boxes; + if (!northPoleTouched) { + // add upper bounding box(es) + if (boxCrosses180Longitude) { + boxes.push_back(box(point(leftBound, upperBound), point(180, 90))); + boxes.push_back(box(point(-180, upperBound), point(rightBound, 90))); + } else { + boxes.push_back(box(point(leftBound, upperBound), point(rightBound, 90))); + } + } + if (!southPoleTouched) { + // add lower bounding box(es) + if (boxCrosses180Longitude) { + boxes.push_back(box(point(leftBound, -90), point(180, lowerBound))); + boxes.push_back(box(point(-180, -90), point(rightBound, lowerBound))); + } else { + boxes.push_back( + box(point(leftBound, -90), point(rightBound, lowerBound))); + } + } + // add the box(es) inbetween the longitude lines + if (boxCrosses180Longitude) { + // only one box needed to cover the longitudes + boxes.push_back(box(point(rightBound, -90), point(leftBound, 90))); + } else { + // two boxes needed, one left and one right of the anti bounding box + boxes.push_back(box(point(-180, -90), point(leftBound, 90))); + boxes.push_back(box(point(rightBound, -90), point(180, 90))); + } + return boxes; +} + +// ____________________________________________________________________________ +bool SpatialJoinAlgorithms::containedInBoundingBoxes( + const std::vector& bbox, point point1) { + // correct lon bounds if necessary + while (point1.get<0>() < -180) { + point1.set<0>(point1.get<0>() + 360); + } + while (point1.get<0>() > 180) { + point1.set<0>(point1.get<0>() - 360); + } + if (point1.get<1>() < -90) { + point1.set<1>(-90); + } else if (point1.get<1>() > 90) { + point1.set<1>(90); + } + + for (size_t i = 0; i < bbox.size(); i++) { + if (boost::geometry::covered_by(point1, bbox.at(i))) { + return true; + } + } + return false; +} + +// ____________________________________________________________________________ +Result SpatialJoinAlgorithms::BoundingBoxAlgorithm() { + const auto [idTableLeft, resultLeft, idTableRight, resultRight, leftJoinCol, + rightJoinCol, numColumns, maxDist, maxResults] = params_; + IdTable result{numColumns, qec_->getAllocator()}; + + // create r-tree for smaller result table + auto smallerResult = idTableLeft; + auto otherResult = idTableRight; + bool leftResSmaller = true; + /*auto smallerChild = childLeft_; + auto otherChild = childRight_; + auto smallerVariable = leftChildVariable_; + auto otherVariable = rightChildVariable_;*/ + auto smallerResJoinCol = leftJoinCol; + auto otherResJoinCol = rightJoinCol; + if (idTableLeft->numRows() > idTableRight->numRows()) { + smallerResult = idTableRight; + otherResult = idTableLeft; + leftResSmaller = false; + /*smallerChild = childRight_; + otherChild = childLeft_; + smallerVariable = rightChildVariable_; + otherVariable = leftChildVariable_; */ + smallerResJoinCol = rightJoinCol; + otherResJoinCol = leftJoinCol; + } + + // Todo in the benchmark: use different algorithms and compare their + // performance + bgi::rtree> rtree; + for (size_t i = 0; i < smallerResult->numRows(); i++) { + // get point of row i + // ColumnIndex smallerJoinCol = getJoinCol(smallerChild, smallerVariable); + auto geopoint = getPoint(smallerResult, i, smallerResJoinCol); + point p(geopoint.value().getLng(), geopoint.value().getLat()); + + // add every point together with the row number into the rtree + rtree.insert(std::make_pair(p, i)); + } + for (size_t i = 0; i < otherResult->numRows(); i++) { + // ColumnIndex otherJoinCol = getJoinCol(otherChild, otherVariable); + auto geopoint1 = getPoint(otherResult, i, otherResJoinCol); + point p(geopoint1.value().getLng(), geopoint1.value().getLat()); + + // query the other rtree for every point using the following bounding box + std::vector bbox = computeBoundingBox(p); + std::vector results; + for (size_t k = 0; k < bbox.size(); k++) { + rtree.query(bgi::intersects(bbox.at(k)), std::back_inserter(results)); + } + for (size_t k = 0; k < results.size(); k++) { + size_t rowLeft = results.at(k).second; + size_t rowRight = i; + if (!leftResSmaller) { + rowLeft = i; + rowRight = results.at(k).second; + } + auto distance = computeDist(idTableLeft, idTableRight, rowLeft, rowRight, + leftJoinCol, rightJoinCol); + if (distance.getDatatype() == Datatype::Int && + distance.getInt() <= maxDist) { + addResultTableEntry(&result, idTableLeft, idTableRight, rowLeft, + rowRight, distance); + } + } + } + Result resTable = + Result(std::move(result), std::vector{}, LocalVocab{}); + return resTable; +} diff --git a/src/engine/SpatialJoinAlgorithms.h b/src/engine/SpatialJoinAlgorithms.h index 00556f1b3..316be773f 100644 --- a/src/engine/SpatialJoinAlgorithms.h +++ b/src/engine/SpatialJoinAlgorithms.h @@ -1,8 +1,20 @@ #pragma once +#include +#include +#include +#include + #include "engine/Result.h" #include "engine/SpatialJoin.h" +namespace bg = boost::geometry; +namespace bgi = boost::geometry::index; + +typedef bg::model::point point; +typedef bg::model::box box; +typedef std::pair value; + class SpatialJoinAlgorithms { public: // initialize the Algorithm with the needed parameters @@ -14,6 +26,16 @@ class SpatialJoinAlgorithms { Result S2geometryAlgorithm(); Result BoundingBoxAlgorithm(); + std::vector OnlyForTestingWrapperComputeBoundingBox( + const point& startPoint) { + return computeBoundingBox(startPoint); + } + + bool OnlyForTestingWrapperContainedInBoundingBoxes( + const std::vector& bbox, point point1) { + return containedInBoundingBoxes(bbox, point1); + } + private: // helper function which returns a GeoPoint if the element of the given table // represents a GeoPoint @@ -33,8 +55,33 @@ class SpatialJoinAlgorithms { const IdTable* resultRight, size_t rowLeft, size_t rowRight, Id distance) const; + // this function computes the bounding box(es), which represent all points, + // which are in reach of the starting point with a distance of at most + // maxDistanceInMeters + std::vector computeBoundingBox(const point& startPoint); + + // this helper function calculates the bounding boxes based on a box, where + // definetly no match can occur. This function gets used, when the usual + // procedure, would just result in taking a big bounding box, which covers + // the whole planet (so for large max distances) + std::vector computeAntiBoundingBox(const point& startPoint); + + // this function returns true, when the given point is contained in any of the + // bounding boxes + bool containedInBoundingBoxes(const std::vector& bbox, point point1); + QueryExecutionContext* qec_; PreparedSpatialJoinParams params_; bool addDistToResult_; std::variant config_; + + // circumference in meters at the equator (max) and the pole (min) (as the + // earth is not exactly a sphere the circumference is different. Note the + // factor of 1000 to convert to meters) + static constexpr double circumferenceMax_ = 40075 * 1000; + static constexpr double circumferenceMin_ = 40007 * 1000; + + // radius of the earth in meters (as the earth is not exactly a sphere the + // radius at the equator has been taken) + static constexpr double radius_ = 6378 * 1000; // * 1000 to convert to meters }; diff --git a/test/engine/SpatialJoinAlgorithmsTest.cpp b/test/engine/SpatialJoinAlgorithmsTest.cpp index f732c0d8e..19382430e 100644 --- a/test/engine/SpatialJoinAlgorithmsTest.cpp +++ b/test/engine/SpatialJoinAlgorithmsTest.cpp @@ -25,7 +25,8 @@ using namespace SpatialJoinTestHelpers; namespace computeResultTest { -class SpatialJoinParamTest : public ::testing::TestWithParam { +class SpatialJoinParamTest + : public ::testing::TestWithParam { public: void createAndTestSpatialJoin( QueryExecutionContext* qec, SparqlTriple spatialJoinTriple, @@ -780,8 +781,277 @@ TEST_P(SpatialJoinParamTest, computeResultSmallDatasetDifferentSizeChildren) { expectedMaxDist10000000_rows_diff, columnNames, false); } -INSTANTIATE_TEST_SUITE_P(SpatialJoin, SpatialJoinParamTest, ::testing::Bool()); +INSTANTIATE_TEST_SUITE_P( + SpatialJoin, SpatialJoinParamTest, + ::testing::Values(SpatialJoin::Algorithm::Baseline, + SpatialJoin::Algorithm::S2Geometry, + SpatialJoin::Algorithm::BoundingBox)); } // end of Namespace computeResultTest +namespace boundingBox { + +namespace bg = boost::geometry; +namespace bgi = boost::geometry::index; + +typedef bg::model::point point; +typedef bg::model::box box; +typedef std::pair value; + +inline void testBoundingBox(const long long& maxDistInMeters, + const point& startPoint) { + auto checkOutside = [&](const point& point1, const point& startPoint, + const std::vector& bbox, + SpatialJoinAlgorithms* spatialJoinAlg) { + // check if the point is contained in any bounding box + bool within = spatialJoinAlg->OnlyForTestingWrapperContainedInBoundingBoxes( + bbox, point1); + if (!within) { + GeoPoint geo1{point1.get<1>(), point1.get<0>()}; + GeoPoint geo2{startPoint.get<1>(), startPoint.get<0>()}; + double dist = ad_utility::detail::wktDistImpl(geo1, geo2) * 1000; + ASSERT_GT(static_cast(dist), maxDistInMeters); + } + }; + + auto testBounds = [](double x, double y, const box& bbox, + bool shouldBeWithin) { + // correct lon bounds if necessary + if (x < -180) { + x += 360; + } else if (x > 180) { + x -= 360; + } + + // testing only possible, if lat bounds are correct and the lon bounds + // don't cover everything (as then left or right of the box is again + // inside the box because of the spherical geometry) + double minLonBox = bbox.min_corner().get<0>(); + double maxLonBox = bbox.max_corner().get<0>(); + if (y < 90 && y > -90 && !(minLonBox < 179.9999 && maxLonBox > 179.9999)) { + bool within = boost::geometry::covered_by(point(x, y), bbox); + ASSERT_TRUE(within == shouldBeWithin); + } + }; + + // build dummy join to access the containedInBoundingBox and + // computeBoundingBox functions. Note that maxDistInMeters has to be accurate, + // otherwise the functions of spatialJoin don't work correctly + auto maxDistInMetersStr = + ""; + auto qec = buildTestQEC(); + auto spatialJoinTriple = + SparqlTriple{TripleComponent{Variable{"?point1"}}, maxDistInMetersStr, + TripleComponent{Variable{"?point2"}}}; + std::shared_ptr spatialJoinOperation = + ad_utility::makeExecutionTree(qec, spatialJoinTriple, + std::nullopt, std::nullopt); + + // add children and test, that multiplicity is a dummy return before all + // children are added + std::shared_ptr op = spatialJoinOperation->getRootOperation(); + SpatialJoin* spatialJoin = static_cast(op.get()); + + PreparedSpatialJoinParams params{ + nullptr, nullptr, nullptr, nullptr, + 0, 0, 1, spatialJoin->getMaxDist(), + std::nullopt}; + + SpatialJoinAlgorithms spatialJoinAlgs{ + qec, params, true, spatialJoin->onlyForTestingGetRealConfig()}; + + std::vector bbox = + spatialJoinAlgs.OnlyForTestingWrapperComputeBoundingBox(startPoint); + // broad grid test + for (int lon = -180; lon < 180; lon += 20) { + for (int lat = -90; lat < 90; lat += 20) { + checkOutside(point(lon, lat), startPoint, bbox, &spatialJoinAlgs); + } + } + + // do tests at the border of the box + for (size_t k = 0; k < bbox.size(); k++) { + // use a small delta for testing because of floating point inaccuracies + const double delta = 0.00000001; + const point minPoint = bbox.at(k).min_corner(); + const point maxPoint = bbox.at(k).max_corner(); + const double lowX = minPoint.get<0>(); + const double lowY = minPoint.get<1>(); + const double highX = maxPoint.get<0>(); + const double highY = maxPoint.get<1>(); + const double xRange = highX - lowX - 2 * delta; + const double yRange = highY - lowY - 2 * delta; + for (size_t i = 0; i <= 100; i++) { + // barely in or out at the left edge + testBounds(lowX + delta, lowY + delta + (yRange / 100) * i, bbox.at(k), + true); + testBounds(lowX - delta, lowY + delta + (yRange / 100) * i, bbox.at(k), + false); + checkOutside(point(lowX - delta, lowY + (yRange / 100) * i), startPoint, + bbox, &spatialJoinAlgs); + // barely in or out at the bottom edge + testBounds(lowX + delta + (xRange / 100) * i, lowY + delta, bbox.at(k), + true); + testBounds(lowX + delta + (xRange / 100) * i, lowY - delta, bbox.at(k), + false); + checkOutside(point(lowX + (xRange / 100) * i, lowY - delta), startPoint, + bbox, &spatialJoinAlgs); + // barely in or out at the right edge + testBounds(highX - delta, lowY + delta + (yRange / 100) * i, bbox.at(k), + true); + testBounds(highX + delta, lowY + delta + (yRange / 100) * i, bbox.at(k), + false); + checkOutside(point(highX + delta, lowY + (yRange / 100) * i), startPoint, + bbox, &spatialJoinAlgs); + // barely in or out at the top edge + testBounds(lowX + delta + (xRange / 100) * i, highY - delta, bbox.at(k), + true); + testBounds(lowX + delta + (xRange / 100) * i, highY + delta, bbox.at(k), + false); + checkOutside(point(lowX + (xRange / 100) * i, highY + delta), startPoint, + bbox, &spatialJoinAlgs); + } + } +} + +TEST(SpatialJoin, computeBoundingBox) { + // ASSERT_EQ("", "uncomment the part below again"); + double circ = 40075 * 1000; // circumference of the earth (at the equator) + // 180.0001 in case 180 is represented internally as 180.0000000001 + for (double lon = -180; lon <= 180.0001; lon += 15) { + // 90.0001 in case 90 is represented internally as 90.000000001 + for (double lat = -90; lat <= 90.0001; lat += 15) { + // circ / 2 means, that all points on earth are within maxDist km of any + // starting point + for (int maxDist = 0; maxDist <= circ / 2.0; maxDist += circ / 36.0) { + testBoundingBox(maxDist, point(lon, lat)); + } + } + } +} + +TEST(SpatialJoin, containedInBoundingBoxes) { + // build dummy join to access the containedInBoundingBox and + // computeBoundingBox functions + auto qec = buildTestQEC(); + auto spatialJoinTriple = SparqlTriple{TripleComponent{Variable{"?point1"}}, + "", + TripleComponent{Variable{"?point2"}}}; + std::shared_ptr spatialJoinOperation = + ad_utility::makeExecutionTree(qec, spatialJoinTriple, + std::nullopt, std::nullopt); + + std::shared_ptr op = spatialJoinOperation->getRootOperation(); + SpatialJoin* spatialJoin = static_cast(op.get()); + + PreparedSpatialJoinParams params{ + nullptr, nullptr, nullptr, nullptr, + 0, 0, 1, spatialJoin->getMaxDist(), + std::nullopt}; + + SpatialJoinAlgorithms spatialJoinAlgs{ + qec, params, true, spatialJoin->onlyForTestingGetRealConfig()}; + + // note that none of the boxes is overlapping, therefore we can check, that + // none of the points which should be contained in one box are contained in + // another box + std::vector boxes = { + box(point(20, 40), point(40, 60)), + box(point(-180, -20), point(-150, 30)), // touching left border + box(point(50, -30), point(180, 10)), // touching right border + box(point(-30, 50), point(10, 90)), // touching north pole + box(point(-45, -90), point(0, -45)) // touching south pole + }; + + // the first entry in this vector is a vector of points, which is contained + // in the first box, the second entry contains points, which are contained in + // the second box and so on + std::vector> containedInBox = { + {point(20, 40), point(40, 40), point(40, 60), point(20, 60), + point(30, 50)}, + {point(-180, -20), point(-150, -20), point(-150, 30), point(-180, 30), + point(-150, 0)}, + {point(50, -30), point(180, -30), point(180, 10), point(50, 10), + point(70, -10)}, + {point(-30, 50), point(10, 50), point(10, 90), point(-30, 90), + point(-20, 60)}, + {point(-45, -90), point(0, -90), point(0, -45), point(-45, -45), + point(-10, -60)}}; + + // all combinations of box is contained in bounding boxes and is not + // contained. a one encodes, the bounding box is contained in the set of + // bounding boxes, a zero encodes, it isn't. If a box is contained, it's + // checked, that the points which should be contained in the box are + // contained. If the box is not contained, it's checked, that the points which + // are contained in that box are not contained in the box set (because the + // boxes don't overlap) + for (size_t a = 0; a <= 1; a++) { + for (size_t b = 0; a <= 1; a++) { + for (size_t c = 0; a <= 1; a++) { + for (size_t d = 0; a <= 1; a++) { + for (size_t e = 0; a <= 1; a++) { + std::vector toTest; // the set of bounding boxes + std::vector> shouldBeContained; + std::vector> shouldNotBeContained; + if (a == 1) { // box nr. 0 is contained in the set of boxes + toTest.push_back(boxes.at(0)); + shouldBeContained.push_back(containedInBox.at(0)); + } else { // box nr. 0 is not contained in the set of boxes + shouldNotBeContained.push_back(containedInBox.at(0)); + } + if (b == 1) { // box nr. 1 is contained in the set of boxes + toTest.push_back(boxes.at(1)); + shouldBeContained.push_back(containedInBox.at(1)); + } else { // box nr. 1 is not contained in the set of boxes + shouldNotBeContained.push_back(containedInBox.at(1)); + } + if (c == 1) { + toTest.push_back(boxes.at(2)); + shouldBeContained.push_back(containedInBox.at(2)); + } else { + shouldNotBeContained.push_back(containedInBox.at(2)); + } + if (d == 1) { + toTest.push_back(boxes.at(3)); + shouldBeContained.push_back(containedInBox.at(3)); + } else { + shouldNotBeContained.push_back(containedInBox.at(3)); + } + if (e == 1) { + toTest.push_back(boxes.at(4)); + shouldBeContained.push_back(containedInBox.at(4)); + } else { + shouldNotBeContained.push_back(containedInBox.at(4)); + } + if (toTest.size() > 0) { + // test all points, which should be contained in the bounding + // boxes + for (size_t i = 0; i < shouldBeContained.size(); i++) { + for (size_t k = 0; k < shouldBeContained.at(i).size(); k++) { + ASSERT_TRUE( + spatialJoinAlgs + .OnlyForTestingWrapperContainedInBoundingBoxes( + toTest, shouldBeContained.at(i).at(k))); + } + } + // test all points, which shouldn't be contained in the bounding + // boxes + for (size_t i = 0; i < shouldNotBeContained.size(); i++) { + for (size_t k = 0; k < shouldNotBeContained.at(i).size(); k++) { + ASSERT_FALSE( + spatialJoinAlgs + .OnlyForTestingWrapperContainedInBoundingBoxes( + toTest, shouldNotBeContained.at(i).at(k))); + } + } + } + } + } + } + } + } +} + +} // namespace boundingBox + } // namespace