Polygon implementation for modern (spherical mean value) interpolation. #329

lewisn-met · 2025-11-05T17:23:55Z

Description

This PR relates to issue #328 and adds:

A SphericalPolygon3D class as an alternative interpolation method to the Triag3D and Quad3D classes currently used in the finite-element interpolation method (FiniteElement.h). This computes weights using the spherical mean value formula rather than the bilinear/barycentric formulae currently implemented.
A test suite for the new class, in test_SphericalPolygon3D.cc, which validates that the interpolation method correctly computes the areas, normals, and weights of an example triangle, planar quadrilateral, and non-planar quadrilateral on the surface of the unit sphere.

(Draft) PR integrating this class into the existing set of interpolation methods in atlas to follow, with the exact implementation to be fixed after further discussion with myself and @odlomax.

Contributor Declaration

By opening this pull request, I affirm the following:

All authors agree to the Contributor License Agreement.
The code follows the project's coding standards.
I have performed self-review and added comments where needed.
I have added or updated tests to verify that my changes are effective and functional.
I have run all existing tests and confirmed they pass.

wdeconinck · 2025-11-07T12:37:08Z

Thanks @lewisn-met for this contribution. I am wondering in how far there is overlap with the ConvexSphericalPolygon which is used in the "ConservativeSphericalPolygon" interpolation, and if missing functionality can be merged perhaps. @sbrdar could you also have a look please?

codecov-commenter · 2025-11-07T13:01:43Z

Codecov Report

❌ Patch coverage is 84.61538% with 32 lines in your changes missing coverage. Please review.
✅ Project coverage is 79.78%. Comparing base (4a5d4e9) to head (70ce3ff).
⚠️ Report is 2 commits behind head on develop.

Files with missing lines	Patch %	Lines
src/tests/interpolation/test_SphericalPolygon3D.cc	80.00%	32 Missing ⚠️

Additional details and impacted files

@@             Coverage Diff             @@
##           develop     #329      +/-   ##
===========================================
+ Coverage    79.27%   79.78%   +0.50%     
===========================================
  Files          909      909              
  Lines        61597    70950    +9353     
===========================================
+ Hits         48832    56605    +7773     
- Misses       12765    14345    +1580

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

lewisn-met · 2025-11-07T14:25:29Z

Thanks for taking a look at this and for the kind words @wdeconinck. From a brief read of the ConvexSphericalPolygon class, the method used to check for the left of a great circle is precisely what I used to check if a point was inside SphericalPolygon3D. As for the intersection method, I will need to look a bit further, but they seem to be distinct.

lewisn-met · 2025-11-07T14:29:39Z

On a slightly unrelated note, I notice that the linux intel-classic ci and gnu-14 ci are are currently failing for this PR. These failures seem to be related to the pluto-benchmark and either the SphericalVector or Cubic2D classes. Is this something you'd like me to take a look at?

wdeconinck · 2025-11-07T15:51:49Z

On a slightly unrelated note, I notice that the linux intel-classic ci and gnu-14 ci are are currently failing for this PR. These failures seem to be related to the pluto-benchmark and either the SphericalVector or Cubic2D classes. Is this something you'd like me to take a look at?

My hunch is that it is due to an updated Eigen version 5.0!!! that happened under the hood without us realising:
linux intel-classic:

  -- Found Eigen3: /home/linuxbrew/.linuxbrew/share/eigen3/cmake/Eigen3Config.cmake (found version "5.0.0")

Any help to fix this is always welcome of course.

wdeconinck · 2025-11-07T20:10:19Z

I have followed up on my hunch, and reverting the CI to use "eigen@3" seems to fix the compilation.
I have added this fix to develop and a 0.44.1 hotfix tag.

This is not a proper fix for compatibility Eigen 5.0.0 but perhaps these are infancy issues that will be sorted upstream in Eigen in due course. At least we can work without failing CI now.
Please rebase on develop.

lewisn-met · 2025-11-10T11:22:02Z

I have now merged the fix from develop into this PR. Thanks for following your hunch so promptly - fingers crossed that any issues with Eigen 5. get sorted out soon!

I also noticed that the validate method included in your ConvexSphericalPolygon class nicely verifies the polygon is on the sphere and points are not too close. I could not see anywhere in the FiniteElement class that does a similar check, so it may be worth sharing this functionality between the methods.

sbrdar · 2025-11-18T11:58:14Z

Hi @lewisn-met. Many thanks for this interesting development. It seems to me that SphericalPolygon3D can be class-derived from ConvexSphericalPolygon. Can we assume that SphericalPolygon3D has to be convex spherical polygon?

lewisn-met · 2025-11-18T14:01:59Z

Hi @sbrdar, thank you for taking a look at this. Indeed, we can assume SphericalPolygon3D is convex, though it need not be (the weights are still positive for star-shaped polygons and well-defined but negative otherwise).

With regards to merging functionality, incorporating the methods of SphericalPolygon3D into ConvexSphericalPolygon, either directly or via inheritance, is certainly an option. The main differences/obstacles here are the representation of the polygon (SphericalPolygon3D requires the triple products from your inLeftHemisphere() function to be stored in addition to the vertices) and encapsulating what is actually necessary for each use case.

Hopefully we can discuss this soon.

sbrdar · 2025-11-19T15:23:42Z

Hi @lewisn-met I am thinking that the Floater's barycentre could be used for the 2nd order cons. remapping in Atlas, so I cherry-picked your implementation with the unit test on SphericalVector3D, and added it to ConvexSphericalPolygon, see the usage here. Maybe something to discuss soon.

lewisn-met · 2025-11-20T11:53:43Z

Would incorporating this method into the cons. remapping be in addition or as an alternative to using Floater's interpolation in the finite-element class?

Our original motivation for the algorithm was as a drop-in alternative to the Quad3D and Triag3D classes used for "continuous" interpolation, see the implementation here. However, I can see how Floater's ideas could be useful in both settings and I don't expect the use of ConvexSphericalPolygon in FiniteElement.cc to be too different from the branch I have linked.

…l-mean-value-polygon

lewisn-met · 2026-01-09T09:05:54Z

@sbrdar I see there is an issue with the preview-publish/deploy action - can you confirm that this is an issue on your end and not with this PR itself? It appears to be related to the python authenticator being called by the job runner.

On a separate note, would you like me to merge from develop, since I see there are several commits which have not yet made their way onto this branch?

sbrdar · 2026-01-09T14:01:33Z

@lewisn-met Thanks very much for the merger with ConvexSphericalPolygon - great work! The tests pass correctly. The failing unrelated static_analyzer problem we can ignore for now and continue with the merge. Can you in the very last step merge from the latest develop?

…cal-mean-value-polygon

lewisn-met · 2026-01-09T14:44:39Z

@sbrdar many thanks for your kind words, I'm glad you like the work! The latest ECMWF/develop has now been merged into this branch and we are good to go.

I have a few small improvements for the spherical mean-value interpolator in mind, but will raise those in a separate PR when I get the chance to work on it. Thanks again for all of your help!

sbrdar

All good for the merge. Very nice results and code implementation. Many thanks!

sbrdar · 2026-01-08T16:42:12Z

src/tests/interpolation/test_interpolation_cubedsphere.cc

+        config.get("type", interpType);
+        config.get("normalisation", normalisationMode);
+
+        SECTION("using " + interpType + " " + normalisationMode) {


A tiny comment here. This output can return, for instance, "using spherical-mean-value true" or "using spherical-mean-value false", but 'false'/'true' refers to normalisation. Can we have normalisationMode replaced with (normalisationMode ? "w/ normalisation" : "")?

sbrdar · 2026-01-08T16:43:45Z

src/tests/interpolation/test_interpolation_cubedsphere.cc

+        config.get("type", interpType);
+        config.get("normalisation", normalisationMode);
+
+        SECTION("using " + interpType + " " + normalisationMode) {


A similar comment as before also here. This output can return, for instance, "using spherical-mean-value true" or "using spherical-mean-value false", but 'false'/'true' refers to normalisation. Can we have normalisationMode replaced with (normalisationMode ? "w/ normalisation" : "")?

sbrdar · 2026-01-08T16:47:09Z

src/tests/interpolation/test_interpolation_cubedsphere.cc

+            }

-    gmshOutput("cubedsphere_source.msh", FieldSet(sourceField));
+            gmshOutput("cubedsphere_source_" + interpType + normalisationMode + ".msh", FieldSet(sourceField));


A tiny suggestion. Perhaps one can replace normalisationMode with (normalisationMode ? "_normalised" : "")

sbrdar · 2026-01-08T16:47:28Z

src/tests/interpolation/test_interpolation_cubedsphere.cc

+            targetFields.add(targetField);
+            targetFields.add(errorField);
+            targetFields.add(partField);
+            gmshOutput("cubedsphere_target_" + interpType + "_" + normalisationMode + ".msh", targetFields);


Here as well. A tiny suggestion. Perhaps one can replace normalisationMode with (normalisationMode ? "_normalised" : "")

sbrdar · 2026-01-08T16:47:48Z

src/tests/interpolation/test_interpolation_finite_element.cc

-    SECTION("test interpolation outputs") {
-        Field field_source = fs.createField<double>(option::name("source"));
-        Field field_target("target", array::make_datatype<double>(), array::make_shape(pointcloud.size()));
+        SECTION("using " + interpType + " " + normalisationMode) {


A tiny suggestion. Perhaps one can replace normalisationMode with (normalisationMode ? "_normalised" : "")

wdeconinck

Hello @lewisn-met I think this is looks a good merge of your previous polygon and this one.
I have some questions and suggestions.

wdeconinck · 2026-01-09T16:43:39Z

src/atlas/util/ConvexSphericalPolygon.cc

 }

+// cf. M. Floater, “Generalized barycentric coordinates and applications” Acta Numerica, p. 001, 2016.
+std::optional<std::vector<double>> ConvexSphericalPolygon::compute_vertex_weights(


Why optional?
Since it is a dynamically sized vector, would it not be same when just passing an empty vector in the "not present" case?
And then in the client code check for vertex_weights.empty()

Actually ignore this, and see below for other suggested API, which I added as a review comment later.

wdeconinck · 2026-01-09T16:54:02Z

src/atlas/util/ConvexSphericalPolygon.cc

+    ATLAS_ASSERT(size_ > 2, "Polygon must have at least 3 points");
+    switch (mode) {
+        case SMV:
+            validateAndComputeNormals();


validateAndComputeNormals() cannot be split into the following, reusing the already available validate?

validate(); compute_edge_normals();

wdeconinck · 2026-01-09T16:55:51Z

src/atlas/util/ConvexSphericalPolygon.h


    template <class Points, typename std::enable_if<contains_PointLonLat<Points>::value, void>::type* = nullptr>
-    ConvexSphericalPolygon(const Points& points): ConvexSphericalPolygon(points.data(), points.size()) {}
+    ConvexSphericalPolygon(const Points& points, const interpolationMode mode = interpolationMode::GRID):


The interpolationMode to me seems fishy.
What is fundamentally different depending on the mode?
The computed area or centroid do not seem different, so why the mode? I would really like to avoid this.

wdeconinck · 2026-01-09T16:57:10Z

src/atlas/util/ConvexSphericalPolygon.h


    size_t size() const { return size_; }

+    std::vector<PointXYZ> edge_normals() const { return edge_normals_; }


This can be computed lazily on demand, just like the area.
Perhaps that could avoid the use of mode.

wdeconinck · 2026-01-09T17:15:25Z

src/atlas/util/ConvexSphericalPolygon.h


+    int previous(const int index) const { return (index == 0) ? size_ - 1 : index - 1; };
+
+    std::optional<std::vector<double>> compute_vertex_weights(const PointXYZ& candidatePoint) const;


Another API may be desired and more performant, avoiding dynamic allocations:

int compute_vertex_weights(const PointXYZ& point, double vertex_weights[], size_t vertex_weights_size) { ATLAS_ASSERT(vertex_weights_size == size()); ... return return_code; // instead of optional }

The intended use would then be compatible with many data container types where the control of allocation is then on the client side.

std::vector<double> vertex_weights; vertex_weights.reserve(max_size); // allocate outside of a loop for (auto polygon : polygons ) { vertex_weights.resize( polygon.size() ); // no reallocation happening if (polygon.compute_vertex_weights( point, vertex_weights.data(), vertex_weights.size() ) == OK){ // success, use vertex_weights } }

or example with std::span (C++20):

std::array<double,1000> buffer; for (auto polygon : polygons ) { std::span<double> vertex_weights(buffer.data(), polygon.size()); if (polygon.compute_vertex_weights( point, vertex_weights.data(), vertex_weights.size() ) == OK) { // success, use vertex_weights } }

lewisn-met added 4 commits November 5, 2025 11:11

Initial commit: added spherical polygon

c119938

Fixed linux build error

aa92414

Debugging test_SphericalPolygon3D

aa78a6f

Debugging rotational consistency test

70ce3ff

github-actions bot added the contributor label Nov 5, 2025

wdeconinck requested a review from sbrdar November 7, 2025 12:34

Merge branch 'ecmwf:develop' into feature/spherical-mean-value-polygon

3d8f72f

lewisn-met added 7 commits November 27, 2025 10:24

Merge remote-tracking branch 'upstream/develop' into feature/spherica…

b058b1a

…l-mean-value-polygon

Merge remote-tracking branch 'upstream/develop' into feature/spherica…

7688593

…l-mean-value-polygon

Added smv option to convexSphericalPolygon constructor.

ca40535

Added interpolation method to ConvexSphericalPolygon

bb7011f

Added weight computation to CSP and tested.

ba58cb0

Fix CI failure - unused variable

8716b30

Added SMV interpolation method and tests

93640db

lewisn-met marked this pull request as draft December 12, 2025 16:08

lewisn-met marked this pull request as ready for review January 6, 2026 10:25

Merge remote-tracking branch 'ECMWFatlas/develop' into feature/spheri…

fad44c7

…cal-mean-value-polygon

sbrdar approved these changes Jan 9, 2026

View reviewed changes

wdeconinck requested changes Jan 9, 2026

View reviewed changes


		size_t size() const { return size_; }

		std::vector<PointXYZ> edge_normals() const { return edge_normals_; }


		int previous(const int index) const { return (index == 0) ? size_ - 1 : index - 1; };

		std::optional<std::vector<double>> compute_vertex_weights(const PointXYZ& candidatePoint) const;

Polygon implementation for modern (spherical mean value) interpolation. #329

Are you sure you want to change the base?

Polygon implementation for modern (spherical mean value) interpolation. #329

Conversation

lewisn-met commented Nov 5, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Contributor Declaration

Uh oh!

wdeconinck commented Nov 7, 2025

Uh oh!

codecov-commenter commented Nov 7, 2025

Codecov Report

Uh oh!

lewisn-met commented Nov 7, 2025

Uh oh!

lewisn-met commented Nov 7, 2025

Uh oh!

wdeconinck commented Nov 7, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

wdeconinck commented Nov 7, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lewisn-met commented Nov 10, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

sbrdar commented Nov 18, 2025

Uh oh!

lewisn-met commented Nov 18, 2025

Uh oh!

sbrdar commented Nov 19, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

lewisn-met commented Nov 20, 2025

Uh oh!

lewisn-met commented Jan 9, 2026

Uh oh!

sbrdar commented Jan 9, 2026

Uh oh!

lewisn-met commented Jan 9, 2026

Uh oh!

sbrdar left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wdeconinck left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

lewisn-met commented Nov 5, 2025 •

edited

Loading

wdeconinck commented Nov 7, 2025 •

edited

Loading

wdeconinck commented Nov 7, 2025 •

edited

Loading

lewisn-met commented Nov 10, 2025 •

edited

Loading

sbrdar commented Nov 19, 2025 •

edited

Loading