-
Notifications
You must be signed in to change notification settings - Fork 43
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
feature: bilinear-remapping interpolation (#85)
* feature: bilinear-remapping interpolation
* fix: make Vector2D constructors public
* fix: missing Vector2D public data method
* fix: public Vector2D x() and y()
* Update src/atlas/interpolation/element/Quad2D.cc
* Update src/atlas/interpolation/element/Quad2D.cc
* Update src/atlas/interpolation/element/Quad2D.cc
* Update src/atlas/interpolation/element/Quad2D.cc
* Update src/atlas/interpolation/element/Quad2D.cc
* Update src/tests/interpolation/test_Triag2D.cc
* Update src/tests/interpolation/test_Triag2D.cc
* Update src/tests/interpolation/test_Triag2D.cc
* Update src/tests/interpolation/test_Triag2D.cc
* Update src/tests/interpolation/test_interpolation_bilinear_remapping.cc
* fix constexpr in Triag2D test
* Apply some more suggestions from code review
* refactor: remove define, use EXPECT_APPROX_EQ
* Fix double definition of BAD_WEIGHT_VALUE
* refactor: response to review
* move logic for quadratic solver into calculation of interpolation
weight
* remove roots struct and unused varibles
* refactor inQuadrilateral
* remove return inside else
* refactor in response to review- Loading branch information
Showing
17 changed files
with
1,813 additions
and
3 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,190 @@ | ||
| /* | ||
| * (C) Crown Copyright 2021 Met Office | ||
| * | ||
| * This software is licensed under the terms of the Apache Licence Version 2.0 | ||
| * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. | ||
| */ | ||
|
|
||
| #include <cmath> | ||
| #include <iostream> | ||
|
|
||
| #include "atlas/interpolation/Vector2D.h" | ||
| #include "atlas/interpolation/element/Quad2D.h" | ||
| #include "atlas/interpolation/element/Triag2D.h" | ||
| #include "atlas/runtime/Log.h" | ||
|
|
||
| namespace atlas { | ||
| namespace interpolation { | ||
| namespace element { | ||
|
|
||
| //---------------------------------------------------------------------------------------------------------------------- | ||
|
|
||
| method::Intersect Quad2D::intersects(const PointXY& r, double edgeEpsilon, double epsilon) const { | ||
| method::Intersect isect; // intersection is false | ||
|
|
||
| /* Split quadrilateral into two triangles, points are labelled counter-clockwise. | ||
| * v01-----------v11 | ||
| * / " . `. | ||
| * / " . `. | ||
| * / " .`. | ||
| * v00-------------------v10 | ||
| * | ||
| */ | ||
|
|
||
| Triag2D T013(v00, v10, v01); | ||
| isect = T013.intersects(r, edgeEpsilon, epsilon); | ||
| if (isect) { | ||
| return isect; | ||
| } | ||
|
|
||
| Triag2D T231(v11, v01, v10); | ||
| isect = T231.intersects(r, edgeEpsilon, epsilon); | ||
| if (isect) { | ||
| isect.u = 1 - isect.u; | ||
| isect.v = 1 - isect.v; | ||
| return isect; | ||
| } | ||
|
|
||
| return isect.fail(); | ||
| } | ||
|
|
||
| method::Intersect Quad2D::localRemap(const PointXY& p, double edgeEpsilon, double epsilon) const { | ||
| method::Intersect isect; | ||
|
|
||
| // work out if point is within the polygon | ||
| if (!inQuadrilateral({p.x(), p.y()})) | ||
| return isect.fail(); | ||
|
|
||
| auto solve_weight = [epsilon](const double a, const double b, const double c, double& wght) -> bool { | ||
| if (std::abs(a) > epsilon) { | ||
| // if a is non-zero, we need to solve a quadratic equation for the weight | ||
| // ax**2 + bx + x = 0 | ||
| double det = b * b - 4. * a * c; | ||
| if (det >= 0.) { | ||
| double inv_two_a = 1. / (2. * a); | ||
| double sqrt_det = std::sqrt(det); | ||
| double root_a = (-b + sqrt_det) * inv_two_a; | ||
| if ((root_a > -epsilon) && (root_a < (1. + epsilon))) { | ||
| wght = root_a; | ||
| return true; | ||
| } | ||
| double root_b = (-b - sqrt_det) * inv_two_a; | ||
| if ((root_b > -epsilon) && (root_b < (1. + epsilon))) { | ||
| wght = root_b; | ||
| return true; | ||
| } | ||
| } | ||
| return false; | ||
| } | ||
| if (std::abs(b) > epsilon) { | ||
| // solve ax + b = 0 | ||
| wght = -c / b; | ||
| return true; | ||
| } | ||
| return false; | ||
| }; | ||
|
|
||
| // solve for u and v where: | ||
| // w1 = ( 1 - u ) * ( 1 - v ) | ||
| // w2 = u * ( 1 - v ) | ||
| // w3 = u * v | ||
| // w4 = ( 1 - u ) * v | ||
|
|
||
| Vector2D ray(p.x(), p.y()); | ||
| Vector2D vA = v00 - ray; | ||
| Vector2D vB = v10 - v00; | ||
| Vector2D vC = v01 - v00; | ||
| Vector2D vD = v00 - v10 - v01 + v11; | ||
|
|
||
| // solve for v | ||
| double a = cross2d(vC, vD); | ||
| double b = cross2d(vC, vB) + cross2d(vA, vD); | ||
| double c = cross2d(vA, vB); | ||
|
|
||
| if (!solve_weight(a, b, c, isect.v)) | ||
| return isect.fail(); | ||
|
|
||
| // solve for u | ||
| a = cross2d(vB, vD); | ||
| b = cross2d(vB, vC) + cross2d(vA, vD); | ||
| c = cross2d(vA, vC); | ||
|
|
||
| if (!solve_weight(a, b, c, isect.u)) | ||
| return isect.fail(); | ||
|
|
||
| return isect.success(); | ||
| } | ||
|
|
||
| bool Quad2D::validate() const { | ||
| // normal for sub-triangle T231 | ||
|
|
||
| Vector2D E23 = v01 - v11; | ||
| Vector2D E21 = v10 - v11; | ||
|
|
||
| double N231 = cross2d(E23, E21); | ||
|
|
||
| // normal for sub-triangle T013 | ||
|
|
||
| Vector2D E01 = v10 - v00; | ||
| Vector2D E03 = v01 - v00; | ||
|
|
||
| double N013 = cross2d(E01, E03); | ||
|
|
||
| // normal for sub-triangle T120 | ||
|
|
||
| Vector2D E12 = -E21; | ||
| Vector2D E10 = -E01; | ||
|
|
||
| double N120 = cross2d(E12, E10); | ||
|
|
||
| // normal for sub-triangle T302 | ||
|
|
||
| Vector2D E30 = -E03; | ||
| Vector2D E32 = -E23; | ||
|
|
||
| double N302 = cross2d(E30, E32); | ||
|
|
||
| // all normals must point same way | ||
|
|
||
| double dot02 = N231 * N013; | ||
| double dot23 = N013 * N120; | ||
| double dot31 = N120 * N302; | ||
| double dot10 = N302 * N231; | ||
|
|
||
| // all normals must point same way | ||
| bool is_inside = ((dot02 >= 0. && dot23 >= 0. && dot31 >= 0. && dot10 >= 0.) || | ||
| (dot02 <= 0. && dot23 <= 0. && dot31 <= 0. && dot10 <= 0.)); | ||
| return is_inside; | ||
| } | ||
|
|
||
| double Quad2D::area() const { | ||
| return std::abs(0.5 * cross2d((v01 - v00), (v11 - v00))) + std::abs(0.5 * cross2d((v10 - v11), (v01 - v11))); | ||
| } | ||
|
|
||
| bool Quad2D::inQuadrilateral(const Vector2D& p) const { | ||
| // point p must be on the inside of all quad edges to be inside the quad. | ||
| return cross2d(p - v00, p - v10) >= 0. && | ||
| cross2d(p - v10, p - v11) >= 0. && | ||
| cross2d(p - v11, p - v01) >= 0. && | ||
| cross2d(p - v01, p - v00) >= 0; | ||
| } | ||
|
|
||
| void Quad2D::print(std::ostream& s) const { | ||
| auto printVector2D = [&s](const Vector2D& v) { s << "[" << v[0] << "," << v[1] << "]"; }; | ||
| s << "Quad2D["; | ||
| printVector2D(v00); | ||
| s << ", "; | ||
| printVector2D(v10); | ||
| s << ", "; | ||
| printVector2D(v11); | ||
| s << ", "; | ||
| printVector2D(v01); | ||
| s << "]"; | ||
| } | ||
|
|
||
|
|
||
| //---------------------------------------------------------------------------------------------------------------------- | ||
|
|
||
| } // namespace element | ||
| } // namespace interpolation | ||
| } // namespace atlas |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,69 @@ | ||
| /* | ||
| * (C) Crown Copyright 2021 Met Office | ||
| * | ||
| * This software is licensed under the terms of the Apache Licence Version 2.0 | ||
| * which can be obtained at http://www.apache.org/licenses/LICENSE-2.0. | ||
| */ | ||
|
|
||
| #pragma once | ||
|
|
||
| #include <iosfwd> | ||
| #include <limits> | ||
|
|
||
| #include "atlas/interpolation/Vector2D.h" | ||
| #include "atlas/interpolation/method/Intersect.h" | ||
| #include "atlas/util/Point.h" | ||
|
|
||
| namespace atlas { | ||
| namespace interpolation { | ||
| namespace method { | ||
| struct Ray; | ||
| } | ||
| namespace element { | ||
|
|
||
| //---------------------------------------------------------------------------------------------------------------------- | ||
|
|
||
| class Quad2D { | ||
| public: | ||
| Quad2D(const double* x0, const double* x1, const double* x2, const double* x3): | ||
| v00(x0), v10(x1), v11(x2), v01(x3) {} | ||
|
|
||
| Quad2D(const PointXY& x0, const PointXY& x1, const PointXY& x2, const PointXY& x3): | ||
| Quad2D(x0.data(), x1.data(), x2.data(), x3.data()) {} | ||
|
|
||
| Quad2D(const Vector2D& x0, const Vector2D& x1, const Vector2D& x2, const Vector2D& x3): | ||
| Quad2D(x0.data(), x1.data(), x2.data(), x3.data()) {} | ||
|
|
||
| method::Intersect intersects(const PointXY& r, double edgeEpsilon = 5 * std::numeric_limits<double>::epsilon(), | ||
| double epsilon = 5 * std::numeric_limits<double>::epsilon()) const; | ||
|
|
||
| method::Intersect localRemap(const PointXY& r, double edgeEpsilon = 5 * std::numeric_limits<double>::epsilon(), | ||
| double epsilon = 5 * std::numeric_limits<double>::epsilon()) const; | ||
|
|
||
| bool validate() const; | ||
|
|
||
| double area() const; | ||
|
|
||
| void print(std::ostream&) const; | ||
|
|
||
| friend std::ostream& operator<<(std::ostream& s, const Quad2D& p) { | ||
| p.print(s); | ||
| return s; | ||
| } | ||
|
|
||
| private: // members | ||
| Vector2D v00; // aka v0 | ||
| Vector2D v10; // aka v1 | ||
| Vector2D v11; // aka v2 | ||
| Vector2D v01; // aka v3 | ||
|
|
||
| static double cross2d(const Vector2D& a, const Vector2D& b) { return a.x() * b.y() - a.y() * b.x(); } | ||
|
|
||
| bool inQuadrilateral(const Vector2D& p) const; | ||
| }; | ||
|
|
||
| //---------------------------------------------------------------------------------------------------------------------- | ||
|
|
||
| } // namespace element | ||
| } // namespace interpolation | ||
| } // namespace atlas |
Oops, something went wrong.