From a45b6e4e6855460b5b2ef4b3b5a1d8924f4c49c5 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Wed, 22 Jan 2025 16:27:39 -0800 Subject: [PATCH 1/5] Contains --- .gitignore | 3 ++- src/base/GridElements.cpp | 12 ++++++++++-- 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/.gitignore b/.gitignore index 315fa2ea..8a54f4b2 100644 --- a/.gitignore +++ b/.gitignore @@ -37,4 +37,5 @@ test/*.eps test/*.png -test/statfigs/* \ No newline at end of file +test/statfigs/* +/.idea diff --git a/src/base/GridElements.cpp b/src/base/GridElements.cpp index bbe8c653..a7dce67a 100644 --- a/src/base/GridElements.cpp +++ b/src/base/GridElements.cpp @@ -144,14 +144,20 @@ bool Face::Contains( const Node & n1 = nodevec[(*this)[i1]]; const Node & n2 = nodevec[(*this)[i2]]; + // TODO: Check if the face contains pole point + + // TODO: Fix the logic here // If both nodes are on the same size of n0.z then there will be no // intersection with the plane z=n0.z. If nodes are on opposite sides // of this plane then they must have an intersection. + + int intersect_flag = 0; // The flag indicating that the edge is very likely to intersect the n0.z plane if ((n1.z > n0.z) && (n2.z > n0.z)) { continue; - } - if ((n1.z < n0.z) && (n2.z < n0.z)) { + }else if ((n1.z < n0.z) && (n2.z < n0.z)) { continue; + } else { + intersect_flag = 1; } // Arcs of constant z aren't informative for determining inside/outside @@ -163,6 +169,8 @@ bool Face::Contains( // Branch here to ensure result is the same regardless of n1-n2 ordering // Under the rules of floating point arithmetic, dA should always be // in the range [0,1]. + + //TODO: incoporate the relative error tolerance from the s2geometry Node nx; if (n1.z < n2.z) { double dA = (n0.z - n1.z) / (n2.z - n1.z); From dd16f3ca4710b0f60730354091e4a37ca0c1241b Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Wed, 22 Jan 2025 20:42:57 -0800 Subject: [PATCH 2/5] remove the "quick check" --- .gitignore | 1 + src/base/GridElements.cpp | 86 +++++++++++++++++++++++++++++++++------ 2 files changed, 74 insertions(+), 13 deletions(-) diff --git a/.gitignore b/.gitignore index 8a54f4b2..96125e6b 100644 --- a/.gitignore +++ b/.gitignore @@ -39,3 +39,4 @@ test/*.png test/statfigs/* /.idea +.vscode/settings.json diff --git a/src/base/GridElements.cpp b/src/base/GridElements.cpp index a7dce67a..ade222da 100644 --- a/src/base/GridElements.cpp +++ b/src/base/GridElements.cpp @@ -138,27 +138,87 @@ bool Face::Contains( ) const { int nParity = 0; + + // TODO: Check if the face contains pole point: Obtain the maximum and the minimum longnitude and check if their difference is larger than Pie + + // Initialize min and max longitude + double minLon = M_PI; + double maxLon = -M_PI; + bool containsPole = false; // Flag to check if the pole might be contained + + // Iterate through edges to calculate min and max longitude + for (size_t i1 = 0; i1 < edges.size(); i1++) { + size_t i2 = (i1 + 1) % edges.size(); + + const Node &n1 = nodevec[edges[i1]]; + const Node &n2 = nodevec[edges[i2]]; + + // Calculate longitude for n1 + double lon1 = atan2(n1.y, n1.x); // Longitude in the range [-π, π] + if (lon1 < 0.0) lon1 += 2 * M_PI; // Convert to [0, 2π] + + // Calculate longitude for n2 + double lon2 = atan2(n2.y, n2.x); // Longitude in the range [-π, π] + if (lon2 < 0.0) lon2 += 2 * M_PI; // Convert to [0, 2π] + + // Update min and max longitude + minLon = std::min({minLon, lon1, lon2}); + maxLon = std::max({maxLon, lon1, lon2}); + } + + // Calculate longitude range + double lonRange = maxLon - minLon; + + + if (lonRange >= M_PI) { + containsPole = true; + + //Now a quick check about the query node location + // If the n0 is the pole point, and the face also contains the same pole point. + if (n0.z == 1 && nodevec[edges[0]][2] < 0.0) { + if (nodevec[edges[0]][2] < 0.0) { + return false; + } else { + return true; + } + } else if (n0.z == -1 && nodevec[edges[0]][2] > 0.0) { + if (nodevec[edges[0]][2] > 0.0) { + return false; + } else { + return true; + } + } + } + + for (size_t i1 = 0; i1 < edges.size(); i1++) { size_t i2 = (i1 + 1) % edges.size(); const Node & n1 = nodevec[(*this)[i1]]; const Node & n2 = nodevec[(*this)[i2]]; - // TODO: Check if the face contains pole point - // TODO: Fix the logic here - // If both nodes are on the same size of n0.z then there will be no - // intersection with the plane z=n0.z. If nodes are on opposite sides - // of this plane then they must have an intersection. - int intersect_flag = 0; // The flag indicating that the edge is very likely to intersect the n0.z plane - if ((n1.z > n0.z) && (n2.z > n0.z)) { - continue; - }else if ((n1.z < n0.z) && (n2.z < n0.z)) { - continue; - } else { - intersect_flag = 1; - } + // Removed the "quick check" that attempted to rule out cases based on both nodes + // being on the same side of n0.z. This approach is not robust and cannot reliably rule out any scenarios. + // In the following, maxLat(n1, n2) refers to the actual maximum latitude along the Great Circle Arc (GCA) + // formed by n1 and n2, including cases where the arc "arches up." + + // 1. If n1.z > n0.z && n2.z > n0.z: + // 1.1 Even if both nodes are above the plane (n0.z),at least one intersection can still occur + // if maxLat(n1, n2) >= n0.z >= minLat(n1, n2). + // + // 2. If n1.z > n0.z > n2.z: + // 2.1 If the GCA does not "arc up," there will be exactly one intersection. + // 2.2 If the GCA "arcs up," and maxLat(n1, n2) occurs due to this arc, then maxLat(n1, n2) > n0.z > n2.z > n1.z. + // In this case, there will be two intersection points. + + // Conclusion: + // The original statement: + // "If both nodes are on the same side of n0.z, then there will be no intersection + // with the plane z = n0.z. If nodes are on opposite sides, they must have an intersection." + // is incorrect and cannot reliably rule out any cases. + // Arcs of constant z aren't informative for determining inside/outside if (n1.z == n2.z) { From 3491c08a9f5b5881d8e3b395e900631eca2c2113 Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Thu, 23 Jan 2025 11:56:26 -0800 Subject: [PATCH 3/5] Update pole query --- src/base/GridElements.cpp | 101 ++++++++++++++++++++------------------ 1 file changed, 54 insertions(+), 47 deletions(-) diff --git a/src/base/GridElements.cpp b/src/base/GridElements.cpp index ade222da..a930c95d 100644 --- a/src/base/GridElements.cpp +++ b/src/base/GridElements.cpp @@ -132,63 +132,70 @@ void Face::RemoveZeroEdges() { /////////////////////////////////////////////////////////////////////////////// +const double DOUBLE_EPS = std::numeric_limits::epsilon(); + + +// Function to classify the location of the face +const int ClassifyFaceLocation(const NodeVector &faceNodes) { + bool hasNorth = false; + bool hasSouth = false; + + for (const auto &node : faceNodes) { + if (node.z > 0) { + hasNorth = true; + } else if (node.z < 0) { + hasSouth = true; + } + + // If both are true, the face crosses hemispheres + if (hasNorth && hasSouth) { + return 0; + } + } + + // Determine if entirely in one hemisphere + if (hasNorth) { + return 1; + } + if (hasSouth) { + return -1; + } + + + return 0; +} + + + bool Face::Contains( const Node & n0, const NodeVector & nodevec ) const { int nParity = 0; - - // TODO: Check if the face contains pole point: Obtain the maximum and the minimum longnitude and check if their difference is larger than Pie + // First check if the query point N0 is around pole area - // Initialize min and max longitude - double minLon = M_PI; - double maxLon = -M_PI; - bool containsPole = false; // Flag to check if the pole might be contained - - // Iterate through edges to calculate min and max longitude - for (size_t i1 = 0; i1 < edges.size(); i1++) { - size_t i2 = (i1 + 1) % edges.size(); - - const Node &n1 = nodevec[edges[i1]]; - const Node &n2 = nodevec[edges[i2]]; - - // Calculate longitude for n1 - double lon1 = atan2(n1.y, n1.x); // Longitude in the range [-π, π] - if (lon1 < 0.0) lon1 += 2 * M_PI; // Convert to [0, 2π] - - // Calculate longitude for n2 - double lon2 = atan2(n2.y, n2.x); // Longitude in the range [-π, π] - if (lon2 < 0.0) lon2 += 2 * M_PI; // Convert to [0, 2π] - - // Update min and max longitude - minLon = std::min({minLon, lon1, lon2}); - maxLon = std::max({maxLon, lon1, lon2}); - } - // Calculate longitude range - double lonRange = maxLon - minLon; + + // Check if the query point is the pole point + if (std::abs(std::abs(n0.z) - 1.0) <= DOUBLE_EPS) { + int face_location = ClassifyFaceLocation(nodevec); + + if (face_location * n0.z > 0 ) { + + + } else if (face_location * n0.z < 0 ) { + return false;// Different hemisphere + + } else { + // The face cross the equator + } + + + } - if (lonRange >= M_PI) { - containsPole = true; - //Now a quick check about the query node location - // If the n0 is the pole point, and the face also contains the same pole point. - if (n0.z == 1 && nodevec[edges[0]][2] < 0.0) { - if (nodevec[edges[0]][2] < 0.0) { - return false; - } else { - return true; - } - } else if (n0.z == -1 && nodevec[edges[0]][2] > 0.0) { - if (nodevec[edges[0]][2] > 0.0) { - return false; - } else { - return true; - } - } - } for (size_t i1 = 0; i1 < edges.size(); i1++) { @@ -203,7 +210,7 @@ bool Face::Contains( // being on the same side of n0.z. This approach is not robust and cannot reliably rule out any scenarios. // In the following, maxLat(n1, n2) refers to the actual maximum latitude along the Great Circle Arc (GCA) // formed by n1 and n2, including cases where the arc "arches up." - + // 1. If n1.z > n0.z && n2.z > n0.z: // 1.1 Even if both nodes are above the plane (n0.z),at least one intersection can still occur // if maxLat(n1, n2) >= n0.z >= minLat(n1, n2). From 8b2167b75895ea15fa1f24e70b18155544d604cf Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Thu, 23 Jan 2025 12:20:19 -0800 Subject: [PATCH 4/5] Update function --- src/base/GridElements.cpp | 42 +++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 10 deletions(-) diff --git a/src/base/GridElements.cpp b/src/base/GridElements.cpp index a930c95d..d85db6b6 100644 --- a/src/base/GridElements.cpp +++ b/src/base/GridElements.cpp @@ -173,23 +173,44 @@ bool Face::Contains( ) const { int nParity = 0; - // First check if the query point N0 is around pole area - - - + // First check if the query point N0 is around pole area. The following check base on the assumption that a face will not contains both ppole points // Check if the query point is the pole point if (std::abs(std::abs(n0.z) - 1.0) <= DOUBLE_EPS) { int face_location = ClassifyFaceLocation(nodevec); - if (face_location * n0.z > 0 ) { + if (face_location * n0.z >= 0 ) {// Same hemisphere + + // Check it's intersection with plane [1,0,0] + // TODO: is the Parity test same for the south hemisphere and cross the equator? + for (size_t i1 = 0; i1 < edges.size(); i1++) { + size_t i2 = (i1 + 1) % edges.size(); + + const Node & n1 = nodevec[(*this)[i1]]; + const Node & n2 = nodevec[(*this)[i2]]; + + // Arcs that go from smaller y to larger y have positive parity. + // Arcs that go from larger y to smaller y have negative parity. + if (n1.y < n2.y) { + nParity++; + + } else { + nParity--; + } + + } + + if (nParity > 0) { + return true; + } else { + return false; + } - } else if (face_location * n0.z < 0 ) { - return false;// Different hemisphere } else { - // The face cross the equator - } + return false;// Different hemisphere + + } } @@ -237,7 +258,8 @@ bool Face::Contains( // Under the rules of floating point arithmetic, dA should always be // in the range [0,1]. - //TODO: incoporate the relative error tolerance from the s2geometry + //TODO: Consider incoporate the relative error tolerance from the s2geometry + // https://github.com/google/s2geometry/blob/d36a49afd22a58be27163d3c835b1e86d312e322/src/s2/s2latlng_rect_bounder.cc#L159 Node nx; if (n1.z < n2.z) { double dA = (n0.z - n1.z) / (n2.z - n1.z); From ee839ebab125936c98c4fda6d2e899885816fb2f Mon Sep 17 00:00:00 2001 From: Hongyu Chen Date: Fri, 24 Jan 2025 11:07:08 -0800 Subject: [PATCH 5/5] Finish comment --- src/base/GridElements.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/base/GridElements.cpp b/src/base/GridElements.cpp index d85db6b6..5ba65ddb 100644 --- a/src/base/GridElements.cpp +++ b/src/base/GridElements.cpp @@ -283,7 +283,10 @@ bool Face::Contains( if (da < 0.0) { continue; - } + } + + // TODO: Break the arc n1 and n2 at nx such that it's mono-tone, + // because even though n1.z < n2.z, it might have two intersection points as I mentioned above // Arcs that go from smaller z to larger z have positive parity. // Arcs that go from larger z to smaller z have negative parity.