feat: Implement is_valid algorithm for polyhedral surfaces #1409

vissarion · 2025-05-21T14:24:14Z

This PR introduces support for the is_valid algorithm for polyhedral surface geometries. A polyhedral_surface is a 3D composite geometry formed by a collection of polygonal faces that meet at shared edges. The implementation follows the OGC Simple Features Specification and introduces structural and topological checks for polyhedral surfaces.

The following properties should hold (following OGC standard):

Contiguity: All polygon patches must form a connected surface via shared boundary segments.
Shared Edges: Shared boundaries must be finite LineStrings used in at most two polygons.
Orientation: Shared boundaries must be traversed in opposite directions by adjacent polygons; global orientation must be consistent.

Six new validity failure types has been introduced that reflect the description of the standard.

Below there is a matrix with valid and invalid cases with visualizations. Those cases (and a few more) are included in the unit tests of this PR.

Name	Geometry	Figure
valid_surface	POLYHEDRALSURFACE(((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((1 1 1,0 1 1,0 0 1,1 0 1,1 1 1)),((1 1 1,1 0 1,1 0 0,1 1 0,1 1 1)),((1 1 1,1 1 0,0 1 0,0 1 1,1 1 1)))
valid_surface_non_convex	POLYHEDRALSURFACE(((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((1 1 1,0 1 1,0 0 1,1 0 1,1 1 1)),((1 1 1,1 0 1,1 0 0,1 1 0,1 .5 .5,1 1 1)))
valid_surface_holes	POLYHEDRALSURFACE(((0 0 0,0 1 0,1 1 0,1 0 0,0 0 0),(0.2 0.2 0,0.2 0.5 0,0.5 0.5 0,0.5 0.2 0,0.2 0.2 0)),((0 0 0,0 0 1,0 1 1,0 1 0,0 0 0)),((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((1 1 1,0 1 1,0 0 1,1 0 1,1 1 1)),((1 1 1,1 0 1,1 0 0,1 1 0,1 1 1)),((1 1 1,1 1 0,0 1 0,0 1 1,1 1 1)))
invalid_surface_non_planar_face	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 1 0,0 0 1,0 0 0)),((1 0 0,1 1 0,0 1 0,1 0 0)))
invalid_surface_collinear_points	POLYHEDRALSURFACE(((1 0 0,1 1 0,0 1 0,1 0 0)),((0 0 0,1 0 0,2 0 0,0 0 0)))
invalid_surface_few_points	POLYHEDRALSURFACE(((1 0 0,1 1 0,0 1 0,1 0 0)),((0 0 0,1 0 0,0 0 0)))
invalid_surface_invalid_intersection_boundary_interior	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 0 1,0 0 0)),((0 1 0,0.5 -0.5 0.5,1 1 0,0 1 0)))
invalid_surface_invalid_intersection_boundary_boundary	POLYHEDRALSURFACE(((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0 -1 1,1 -1 1,1 1 0,0 1 0)))
invalid_surface_invalid_intersection_in_common_vertex	POLYHEDRALSURFACE(((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,1 0 1,1 1 0,0 1 0)))
invalid_surface_invalid_intersection_partial_edge	POLYHEDRALSURFACE(((0 0 0,1 0 0,1 0 1,0 0 1,0 0 0)),((0 1 0,0.5 0 1,1 0 1,1 1 0,0 1 0)))
invalid_surface_disconnected_polygons	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 1 0,0 0 0)),((1 1 0,2 1 0,2 2 0,1 2 0,1 1 0)))
invalid_surface_invalid_intersection_vertex_edge	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 1 0,0 0 0)),((0.5 0.5 0,2 1 0,2 2 0,1 2 0,0.5 0.5 0)))
invalid_surface_invalid_intersection_vertex_vertex	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 1 0,0 0 0)),((1 0 0,2 1 0,2 2 0,1 2 0,1 0 0)))
invalid_surface_invalid_intersection_parallel_edges	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 1 0,0 0 0)),((0 1 0,0.5 0.5 0,1 0 0,2 1 0,2 2 0,1 2 0,0 1 0)))
invalid_surface_invalid_intersection_operlapping_faces	POLYHEDRALSURFACE(((0 0 0,1 0.5 0,0.5 1 0,0 0 0)),((0 1 0,1 0 0,2 1 0,2 2 0,1 2 0,0 1 0)))
invalid_surface_invalid_intersection_operlapping_faces2	POLYHEDRALSURFACE(((0 0 0,1 0 0,0.5 1 0,0 0 0)),((0 1 0,1 0 0,2 1 0,2 2 0,1 2 0,0 1 0)))
invalid_surface_invalid_face	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 1 0,0 0 0)),((0 1 0,1 0 0,2 2 0,2 1 0,1 2 0,0 1 0)))
invalid_surface_non_manifold_edge	POLYHEDRALSURFACE(((0 0 0,1 0 0,0 1 0,0 0 0)),((0 1 0,1 0 0,2 1 0,2 2 0,1 2 0,0 1 0)),((0 1 0,1 0 0,1 1 1,0 1 0)))

Two unit tests are failing and are commented out. Those are dependent on #1406

Copilot

Pull Request Overview

This PR enhances the validity checking algorithm by introducing support for polyhedral surfaces in line with the OGC Simple Features Specification. Key changes include:

Adding several new failure types for polyhedral surfaces (e.g. collinear, non-coplanar points, inconsistent orientation, etc.).
Extending the test suite in is_valid.cpp with comprehensive cases for various valid and invalid polyhedral surfaces.
Updating include directives to incorporate polyhedral surface support in the core algorithm implementation.

Reviewed Changes

Copilot reviewed 4 out of 4 changed files in this pull request and generated 2 comments.

File	Description
test/algorithms/is_valid.cpp	Adds polyhedral surface test cases and minor documentation updates.
include/boost/geometry/algorithms/validity_failure_type.hpp	Introduces new failure codes specific to polyhedral surfaces.
include/boost/geometry/algorithms/detail/is_valid/implementation.hpp	Adds header inclusion for polyhedral surface implementation.

test/algorithms/is_valid.cpp

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

tinko92

I noted some points inline, mostly minor but importantly (imho) the question of separating some computations into strategy.

A minor point, this PR should probably include the removal of

geometry/include/boost/geometry/geometries/polyhedral_surface.hpp

Line 47 in d2e5137

//TODO: add is_valid() check

tinko92 · 2025-05-21T15:19:32Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+    static bool are_linearly_dependent(vec3 const& a, vec3 const& b)
+    {
+        return qvm::mag_sqr(qvm::cross(a, b)) < tolerance * tolerance;
+    }
+
+    // Try to find 3 non-collinear points
+    static bool find_non_collinear_points(polygon_3d const& polygon, vec3& base, vec3& v1, vec3& v2)
+    {
+        for (auto pit1 = points_begin(polygon); pit1 != points_end(polygon); ++pit1)
+        {
+            for (auto pit2 = pit1; pit2 != points_end(polygon); ++pit2)
+            {
+                for (auto pit3 = pit2; pit3 != points_end(polygon); ++pit3)
+                {
+                    set_vec_from_point::apply(*pit1, base);
+                    set_vec_from_point::apply(*pit2, v1);
+                    set_vec_from_point::apply(*pit3, v2);
+                    if (!are_linearly_dependent(v1 - base, v2 - base))
+                    {
+                        return true;
+                    }
+                }
+            }
+        }
+        return false; // All points are collinear
+    }


I think by the library design philosophy and also practically, if one wants to be able to plug in a robust computation later, this is a cartesian strategy.

Indeed, a lot of computations here are coordinate system spesific and should be moved in strategies. On one hand, I think that polyhedral surface can only be cartesian. But on the other hand, predicates should be moved to stategies anyways for robustness reasons (as you said).

I thought about it. I think, these are equivalent to three calls to a side_strategy for the the XY, XZ, YZ projections respectively. If any of them returns non-zero, they are non-collinear. We already have the side strategy with the correct tolerances for this,

geometry/include/boost/geometry/strategies/cartesian/side_rounded_input.hpp

Line 34 in d2e5137

struct side_rounded_input

.

This can short-circuit (if the call for the first projection returns non-zero which it will almost surely in general position, then the other calls are not necessary) and there is probably a lot of common subexpression elimination that the compiler will do.

That reduces the computation to orientation, great. I was thinking using the project_point_to_plane_2d to only do one orientation test but I think projecting to coordinate planes as you suggested is better since it does not involve any computations (substitutions, dot products etc).
I change this part to use the side strategy provided by the umbrella strategy used in is_valid. This uses the side_by_triangle strategy as the default but the user can change it (as in many other cases in the library). At this point I think this is better than creating a new umbrella strategy that has side_rounded_input as a side strategy. We can think of substituting the side_by_triangle by side_rounded_input in general or in as many algorithms as possible as a future step.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

tinko92 · 2025-05-21T15:49:33Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+                vec3 vi;
+                set_vec_from_point::apply(*pit, vi);
+                auto vol = qvm::dot(normal, vi)- dot(normal, base);
+                if (std::abs(vol) > tolerance)


This may be to narrow for e.g. float coordinates in the range of thousand, I think, the bound should probably be scaled to the input magnitude, I will revisit this later when I have more time.

This is orient3d (see e.g. https://www.[cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c](https://www.cs.cmu.edu/afs/cs/project/quake/public/code/predicates.c)), i.e. the 3D generalisation of side, a determinant involving *pit, base, v1, v2. Bound for a tolerant version can be derived the same way as for side_rounded_input, I'll try to do it by the end of the week. Calling orient3d for every point may look more expensive, but I think the compiler will know to pull most unchanging subexpressions out of the loop (same as dot(normal, base) could be pulled out of the loop).

Exactly, orient3d is equivalent. I think it would be useful to add some explanations and/or references on how the tolerance bounds are derived as comments in side_rounded_input.

I thought about it. I think, we want
$$\begin{vmatrix} a_x-d_x & a_y -d_y & a_z - d_z \\\ b_x-d_x & b_y-d_y & b_z-d_z \\\ c_x-d_x & c_y-d_y & c_z-d_z \end{vmatrix}$$
Since only one point changes (let's say $a=\text{vi}$), it makes sense to take the first row as cofactors and the minors will be the elements of $n$ (if all points including $a$ were translated by -base before). Then the determinant is

$(a_x-d_x) \times ((b_y-d_y) \times (c_z - d_z) - (b_z - d_z) \times (c_y - d_y)) - (a_y - d_y) \times \ldots + (a_z - d_z) \times \ldots$

Then, I think, the tolerance (based on https://zenodo.org/records/7539355 with rounded input rules, rounding the second order part up because it doesn't matter for this use case):
$$12\epsilon \times \left(\left(\left\lvert a_x\right\rvert + \left\lvert d_x\right\rvert \right) \times \left(\left(\left\lvert b_y\right\rvert + \left\lvert d_y\right\rvert \right) \times \left(\left\lvert c_z\right\rvert + \left\lvert d_z\right\rvert \right) + \left(\left\lvert b_z\right\rvert + \left\lvert d_z\right\rvert \right) \times \left(\left\lvert c_y\right\rvert + \left\lvert d_y\right\rvert \right)\right) + \left(\left\lvert a_y\right\rvert + \left\lvert d_y\right\rvert \right) \times \left(\left(\left\lvert b_x\right\rvert + \left\lvert d_x\right\rvert \right) \times \left(\left\lvert c_z\right\rvert + \left\lvert d_z\right\rvert \right) + \left(\left\lvert b_z\right\rvert + \left\lvert d_z\right\rvert \right) \times \left(\left\lvert c_x\right\rvert + \left\lvert d_x\right\rvert \right)\right) + \left(\left\lvert a_z\right\rvert + \left\lvert d_z\right\rvert \right) \times \left(\left(\left\lvert b_x\right\rvert + \left\lvert d_x\right\rvert \right) \times \left(\left\lvert c_y\right\rvert + \left\lvert d_y\right\rvert \right) + \left(\left\lvert b_y\right\rvert + \left\lvert d_y\right\rvert \right) \times \left(\left\lvert c_x\right\rvert + \left\lvert d_x\right\rvert \right)\right)\right)$$
where $\epsilon$ is half of numerical_limits::epsilon(). This is the same structure as the expanded determinant. If I made no mistake, this will guarantee that the tolerance is sufficiently large such that for all input points that are correctly rounded approximations of actually collinear points, it will find them coplanar.

Looks like a lot of ops but since only a changes, most of the subexpressions should be eliminated in the loop and it should just be around something like ~22 math ops per point.

And this would be the equivalent of side_rounded for planes.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-05-26T16:49:29Z

Looks great! I will review it next Wednesday or Thursday.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

vissarion · 2025-05-28T12:57:44Z

@tinko92 thanks for the comments. I fixed all but the ones about tolerance and the strategy. I need to think of a unified and consistent solution for those.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-05-28T12:58:52Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+            typename boost::geometry::model::polygon<point_2d, false>
+        >::type;
+    using ring_2d = typename boost::geometry::ring_type<polygon_2d>::type;
+    using segment_2d = typename boost::geometry::model::segment<point_2d>;


Do we need this many helper types?

segment_2d is not used any more so I removed. Ring types we can get them from std::dectype if you prefer it, I am ok with those as is.

It's fine. But did you push already? I don't see the updates...

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-05-28T13:15:23Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        }
+        else
+        {
+            auto const& ring = boost::geometry::exterior_ring(polygon);


So this block does really the projection. Nice.
I think it's good to extract it from this class - and even to another source.
Such that it is useful on its own, and can be unit tested on its own.

I agree, would you mind leaving it as future work to see also what other algorithms could use it? As it is now there are a lot of shared variables with the rest of the code as well as errors related to is_valid.

But exactly that (shared variables) will makes it harder to extract it later as well...

A generic polygon projection from 3d to 2d plane would be great. We could then immediately implement the surface-area algorithm on a 3d surface...

It would be a generic algorithm, every patch (surface) would be a separate 2D polygon. Or is it different than what you do here (I couldn't look into all the details yet).

Anyway - I'm fine having it as a next PR - I'm only afraid that it would not really safe time.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-05-28T13:16:43Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+                    return intersection_type::invalid;
+                }
+            }
+            // Project the polygon onto the 2D plane defined by v1 and v2


Looks interesting. I'm not experienced in this so cannot really comment.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-05-28T13:20:03Z

Looks great! I will review it next Wednesday or Thursday.

I will continue tomorrow. It looks good but I had many suggestions already...

tinko92 · 2025-05-28T14:03:27Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+                    }
+                    else if (intersection == intersection_type::valid)
+                    {
+                        disconnected = false;


Does this check work? My understanding from the specification is that the entire polyhedral surface must be connected, is this not different from any pair of polygons being connected? If my understanding of the standard is correct, this needs a BFS.

You are right. Currently, the check only checks that there is no disconnected polygon. More work is needed to check for connectedness of the surface.

barendgehrels · 2025-06-04T18:20:49Z

Looks great! I will review it next Wednesday or Thursday.

I will continue tomorrow. It looks good but I had many suggestions already...

So I didn't continue tomorrow - but great that you processed our comments! Thanks! I plan now to continue this weekend.

vissarion · 2025-06-04T19:31:17Z

I refactor the code addressing almost all of the comments. Some things are left as todo for next PR(s). What is missing from this PR is the predicates (places in the code that use tolerance) that should be written as orientation checks implemented in strategies. I will work on it in the next days.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-06-08T07:41:49Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        return strategy.side().apply(p1_2, p2_2, p3_2) != 0;
+    }
+
+    // Try to find 3 non-collinear points


in any plane?

3 non-colinear points will define a plane. Later in the code we check if the rest of the points of the polygon lay on this plane.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-06-08T07:57:07Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        return false; // All points are collinear
+    }
+
+    static point_2d project_point_to_plane_2d(vec3 const& v1, vec3 const& v2, vec3 const& p)


Without comments, and drawings, and much experience, it's hard to follow this...

barendgehrels · 2025-06-08T07:58:05Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        auto ring = boost::geometry::exterior_ring(polygon);
+        if (boost::geometry::closure<ring_2d>::value == boost::geometry::closed)
+        {
+            ring.pop_back(); // Remove the closing point to make it open


We have also an iterator doing this automatically. It closes it automatically (which is similar - but you can then always skip the last point).
That avoids the need for a copy of the ring.

Thanks, I will use closing_iterator.

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

barendgehrels · 2025-06-08T08:02:38Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        return point_2d{qvm::dot(v1, p), qvm::dot(v2, p)};
+    }
+
+    template <typename VisitPolicy>


So what does it do exactly? It finds if the segment point1-point2 overlaps with any segment of the outer ring of the polygon?

To check intersection of a segment point1-point2 and a polygon (all in 3D) it computes the intersection of the segment with the plane of the the polygon and then checks whether this projected intersection point overlaps with the projected polygon.

barendgehrels · 2025-06-08T08:04:25Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        // Check if the polygon has been projected before
+        if (!boost::geometry::is_empty(m_polygon_projections[polygon_index].projected_polygon))
+        {
+            auto const& projection = m_polygon_projections[polygon_index];


Shouldn't we project all polygons? Are there cases where projecting a polygon is unnecessary work?
If we project all of them, we could do that at the beginning.

But probably there are conditions that we can stop already, so then this approach is better indeed. That could use a comment.

barendgehrels · 2025-06-08T08:05:04Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        else
+        {
+            // TODO: Support interior rings
+            ring_3d const ring = boost::geometry::exterior_ring(polygon);


It's a copy...

barendgehrels · 2025-06-08T08:17:29Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+    using ring_3d = typename boost::geometry::ring_type<polygon_3d>::type;
+    using point_3d = typename boost::geometry::point_type<ring_3d>::type;
+    using segment_3d = typename boost::geometry::model::segment<point_3d>;
+    using set_vec_from_point = typename strategy::transform::detail::matrix_transformer::


The name is not 100% clear to me.

But (I see now) we have that already for a long time.
It sets a vector from a point?
Anyway, the using clause is fine then...

barendgehrels · 2025-06-08T08:18:25Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        return false;
+    }
+
+    // Check if a segment intersects a polygon in 3D space


So also here point1 and point2 form a segment, it would be good to make that clear.

barendgehrels · 2025-06-08T08:19:18Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+    intersection_type check_polygon_segment_intersection(point_3d const& point1,
+                                                         point_3d const& point2,
+                                                         polygon_3d const& polygon,
+                                                         size_t const& polygon_index,


-> std::size_t and for this one the & is not necessary (neither harmful but we never do it). const we often do neither but it would actually be good so that can be kept.

barendgehrels · 2025-06-08T08:20:17Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        vec3 p1;
+        vec3 p2;
+        set_vec_from_point::apply(point1, p1);
+        set_vec_from_point::apply(point2, p2);


Why not v1 and v2 here, if these are vectors?
And then below you have another set of v1 and v2, it's not clear to me...

barendgehrels · 2025-06-08T08:21:44Z

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp

+        // Check if the segment intersects the plane of the polygon
+        coordinate_type const dot_product = dot(normal, base);
+        coordinate_type const sign1 = dot(normal, p1) - dot_product;
+        coordinate_type const sign2 = dot(normal, p2) - dot_product;


The p1 and p2 are calculated long before...
Would it be possible to move that to here? (And rename them to something vector-like)?

vissarion · 2025-06-12T15:22:14Z

@barendgehrels thanks for you comments! I addressed most of them. I need to work more on this PR. Apart from the orient3d some other simplifications can be implemented (the unaddressed comments of Barend are related to those simplifications). I will work on it next week.

vissarion requested review from barendgehrels, awulkiew, tinko92 and Copilot May 21, 2025 14:24

Copilot AI reviewed May 21, 2025

View reviewed changes

test/algorithms/is_valid.cpp Show resolved Hide resolved

test/algorithms/is_valid.cpp Show resolved Hide resolved

tinko92 reviewed May 21, 2025

View reviewed changes

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp Outdated Show resolved Hide resolved

tinko92 reviewed May 21, 2025

View reviewed changes

tinko92 reviewed May 25, 2025

View reviewed changes

include/boost/geometry/algorithms/detail/is_valid/polyhedral_surface.hpp Outdated Show resolved Hide resolved

barendgehrels assigned vissarion May 26, 2025

barendgehrels added the enhancement label May 26, 2025

vissarion force-pushed the feat/polysurf_is_valid branch 2 times, most recently from fe6743a to 36cb224 Compare May 28, 2025 12:35

barendgehrels reviewed May 28, 2025

View reviewed changes