/src/geos/src/operation/grid/TraversalAreas.cpp

Source
/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2018-2025 ISciences, LLC
 *
 * This is free software; you can redistribute and/or modify it under
 * the terms of the GNU Lesser General Public Licence as published
 * by the Free Software Foundation.
 * See the COPYING file for more information.
 *
 **********************************************************************/

#include <limits>
#include <vector>

#include <geos/algorithm/Area.h>
#include <geos/algorithm/Orientation.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/operation/grid/PerimeterDistance.h>

#include <geos/operation/polygonize/Polygonizer.h>
#include <geos/operation/grid/TraversalAreas.h>

using geos::geom::CoordinateSequence;
using geos::geom::CoordinateXY;
using geos::geom::Geometry;
using geos::geom::GeometryFactory;
using geos::geom::Envelope;

namespace geos::operation::grid {

struct CoordinateChain
{
    double start; // perimeter distance value of the first coordinate
    double stop;  // perimeter distance value of the last coordinate
    const std::vector<CoordinateXY>* coordinates;
    bool visited;

    CoordinateChain(double p_start, double p_stop, const std::vector<CoordinateXY>* p_coords)
      : start{ p_start }
      , stop{ p_stop }
      , coordinates{ p_coords }
      , visited{ false }
    {
    }
};

static double
exit_to_entry_perimeter_distance_ccw(const CoordinateChain& c1, const CoordinateChain& c2, double perimeter)
{
    return PerimeterDistance::getPerimeterDistanceCCW(c1.stop, c2.start, perimeter);
}

static CoordinateChain*
getNextChain(std::vector<CoordinateChain>& chains,
             const CoordinateChain* chain,
             const CoordinateChain* kill,
             double perimeter)
{

    CoordinateChain* min = nullptr;
    double min_distance = std::numeric_limits<double>::max();

    for (CoordinateChain& candidate : chains) {
        if (candidate.visited && std::addressof(candidate) != kill) {
            continue;
        }

        double distance = exit_to_entry_perimeter_distance_ccw(*chain, candidate, perimeter);
        if (distance < min_distance) {
            min_distance = distance;
            min = std::addressof(candidate);
        }
    }

    if (min == nullptr) {
        throw std::runtime_error("Failed to find next CoordinateChain");
    }

    return min;
}

template<typename T>
bool
hasMultipleUniqueCoordinates(const T& vec)
{
    for (std::size_t i = 1; i < vec.size(); i++) {
        if (vec[i] != vec[0]) {
            return true;
        }
    }

    return false;
}

/**
 * @brief Identify counter-clockwise rings formed by the supplied coordinate vectors and the boundary of this box.
 *
 * @param box Box to be included in rings
 * @param coord_lists A list of coordinate vectors, with each vector representing a traversal of `box` or a closed ring.
 * @param visitor Function be applied to each ring identified. Because `coord_lists` may include both clockwise and
 *                counter-clockwise closed rings, `visitor` will be provided with the orientation of each ring as an
 *                argument.
 */
template<typename F>
void
visitRings(const Envelope& box, const std::vector<const std::vector<CoordinateXY>*>& coord_lists, F&& visitor)
{
    std::vector<CoordinateChain> chains;
    chains.reserve(coord_lists.size() + 4);

    for (const auto& coords : coord_lists) {
        if (!hasMultipleUniqueCoordinates(*coords)) {
            continue;
        }

        if (coords->front() == coords->back() && hasMultipleUniqueCoordinates(*coords)) {
            // Closed ring. Check orientation.

            // TODO: Remove copy
            CoordinateSequence seq(0, false, false);
            seq.setPoints(*coords);
            bool is_ccw = algorithm::Orientation::isCCW(&seq);
            visitor(*coords, is_ccw);
        } else {
            double start = PerimeterDistance::getPerimeterDistance(box, coords->front());
            double stop = PerimeterDistance::getPerimeterDistance(box, coords->back());

            chains.emplace_back(start, stop, coords);
        }
    }

    double height{ box.getHeight() };
    double width{ box.getWidth() };
    double perimeter{ box.getPerimeter() };

    // create coordinate lists for corners
    std::vector<CoordinateXY> bottom_left = { CoordinateXY(box.getMinX(), box.getMinY()) };
    std::vector<CoordinateXY> top_left = { CoordinateXY(box.getMinX(), box.getMaxY()) };
    std::vector<CoordinateXY> top_right = { CoordinateXY(box.getMaxX(), box.getMaxY()) };
    std::vector<CoordinateXY> bottom_right = { CoordinateXY(box.getMaxX(), box.getMinY()) };

    // Add chains for corners
    chains.emplace_back(0.0, 0.0, &bottom_left);
    chains.emplace_back(height, height, &top_left);
    chains.emplace_back(height + width, height + width, &top_right);
    chains.emplace_back(2 * height + width, 2 * height + width, &bottom_right);

    std::vector<CoordinateXY> coords;
    for (auto& chain_ref : chains) {
        if (chain_ref.visited || chain_ref.coordinates->size() == 1) {
            continue;
        }

        coords.clear();

        CoordinateChain* chain = std::addressof(chain_ref);
        CoordinateChain* first_chain = chain;
        do {
            chain->visited = true;
            coords.insert(coords.end(), chain->coordinates->cbegin(), chain->coordinates->cend());
            chain = getNextChain(chains, chain, first_chain, perimeter);
        } while (chain != first_chain);

        coords.push_back(coords[0]);

        if (hasMultipleUniqueCoordinates(coords)) {
            visitor(coords, true);
        }
    }
}

double
TraversalAreas::getLeftHandArea(const Envelope& box, const std::vector<const std::vector<CoordinateXY>*>& coord_lists)
{
    double ccw_sum = 0;
    double cw_sum = 0;
    bool found_a_ring = false;

    visitRings(box, coord_lists, [&cw_sum, &ccw_sum, &found_a_ring](const std::vector<CoordinateXY>& coords, bool is_ccw) {
        found_a_ring = true;

        if (is_ccw) {
            ccw_sum += algorithm::Area::ofRing(coords);
        } else {
            cw_sum += algorithm::Area::ofRing(coords);
        }
    });

    if (!found_a_ring) {
        throw std::runtime_error("Cannot determine coverage fraction (it is either 0 or 100%)");
    }

    // If this box has only clockwise rings (holes) then the area
    // of those holes should be subtracted from the complete box area.
    if (ccw_sum < cw_sum) {
        ccw_sum += box.getArea();
    }

    return ccw_sum - cw_sum;
}

std::unique_ptr<Geometry>
TraversalAreas::getLeftHandRings(const GeometryFactory& gfact, const Envelope& box, const std::vector<const std::vector<CoordinateXY>*>& coord_lists)
{
    using geom::LinearRing;

    std::vector<std::unique_ptr<LinearRing>> shells;
    std::vector<std::unique_ptr<LinearRing>> holes;

    bool found_a_ring = false;

    visitRings(box, coord_lists, [&gfact, &shells, &holes, &found_a_ring](const std::vector<CoordinateXY>& coords, bool is_ccw) {
        found_a_ring = true;

        // finding a collapsed ring is sufficient to determine whether the cell interior is inside our outside,
        // but we don't want to actually construct the ring.
        if (algorithm::Area::ofRing(coords) == 0) {
            return;
        }

        CoordinateSequence seq(0, false, false);
        seq.setPoints(coords);
        auto ring = gfact.createLinearRing(std::move(seq));

        if (is_ccw) {
            shells.push_back(std::move(ring));
        } else {
            holes.push_back(std::move(ring));
        }
    });

    if (!found_a_ring) {
        throw std::runtime_error("Cannot determine coverage fraction (it is either 0 or 100%)");
    }

    if (shells.empty() && holes.empty()) {
        return gfact.createPolygon();
    }

    if (shells.empty() && !holes.empty()) {
        CoordinateSequence seq(5, false, false);
        seq.setAt(CoordinateXY{box.getMinX(), box.getMinY()}, 0);
        seq.setAt(CoordinateXY{box.getMaxX(), box.getMinY()}, 1);
        seq.setAt(CoordinateXY{box.getMaxX(), box.getMaxY()}, 2);
        seq.setAt(CoordinateXY{box.getMinX(), box.getMaxY()}, 3);
        seq.setAt(CoordinateXY{box.getMinX(), box.getMinY()}, 4);

        shells.push_back(gfact.createLinearRing(std::move(seq)));
    }

    if (holes.empty()) {
        std::vector<std::unique_ptr<Geometry>> polygons;
        for (auto& shell : shells) {
            polygons.push_back(gfact.createPolygon(std::move(shell)));
        }

        if (polygons.size() == 1) {
            return std::move(polygons.front());
        } else {
            return gfact.createMultiPolygon(std::move(polygons));
        }
    }

    if (shells.size() == 1) {
        return gfact.createPolygon(std::move(shells.front()), std::move(holes));
    }

    polygonize::Polygonizer polygonizer(true);

    for (const auto& shell : shells) {
      polygonizer.add(static_cast<const Geometry*>(shell.get()));
    }
    for (const auto& hole : holes) {
      polygonizer.add(static_cast<const Geometry*>(hole.get()));
    }

    auto polygonized = polygonizer.getPolygons();
    return gfact.createMultiPolygon(std::move(polygonized));
}

}

Coverage Report

Created: 2025-10-28 07:10

Line	Count	Source
1		/**********************************************************************
2		*
3		* GEOS - Geometry Engine Open Source
4		* http://geos.osgeo.org
5		*
6		* Copyright (C) 2018-2025 ISciences, LLC
7		*
8		* This is free software; you can redistribute and/or modify it under
9		* the terms of the GNU Lesser General Public Licence as published
10		* by the Free Software Foundation.
11		* See the COPYING file for more information.
12		*
13		**********************************************************************/
14
15		#include <limits>
16		#include <vector>
17
18		#include <geos/algorithm/Area.h>
19		#include <geos/algorithm/Orientation.h>
20		#include <geos/geom/Coordinate.h>
21		#include <geos/geom/CoordinateSequence.h>
22		#include <geos/geom/Geometry.h>
23		#include <geos/geom/GeometryFactory.h>
24		#include <geos/operation/grid/PerimeterDistance.h>
25
26		#include <geos/operation/polygonize/Polygonizer.h>
27		#include <geos/operation/grid/TraversalAreas.h>
28
29		using geos::geom::CoordinateSequence;
30		using geos::geom::CoordinateXY;
31		using geos::geom::Geometry;
32		using geos::geom::GeometryFactory;
33		using geos::geom::Envelope;
34
35		namespace geos::operation::grid {
36
37		struct CoordinateChain
38		{
39		double start; // perimeter distance value of the first coordinate
40		double stop; // perimeter distance value of the last coordinate
41		const std::vector<CoordinateXY>* coordinates;
42		bool visited;
43
44		CoordinateChain(double p_start, double p_stop, const std::vector<CoordinateXY>* p_coords)
45	0	: start{ p_start }
46	0	, stop{ p_stop }
47	0	, coordinates{ p_coords }
48	0	, visited{ false }
49	0	{
50	0	}
51		};
52
53		static double
54		exit_to_entry_perimeter_distance_ccw(const CoordinateChain& c1, const CoordinateChain& c2, double perimeter)
55	0	{
56	0	return PerimeterDistance::getPerimeterDistanceCCW(c1.stop, c2.start, perimeter);
57	0	}
58
59		static CoordinateChain*
60		getNextChain(std::vector<CoordinateChain>& chains,
61		const CoordinateChain* chain,
62		const CoordinateChain* kill,
63		double perimeter)
64	0	{
65
66	0	CoordinateChain* min = nullptr;
67	0	double min_distance = std::numeric_limits<double>::max();
68
69	0	for (CoordinateChain& candidate : chains) {
70	0	if (candidate.visited && std::addressof(candidate) != kill) {
71	0	continue;
72	0	}
73
74	0	double distance = exit_to_entry_perimeter_distance_ccw(*chain, candidate, perimeter);
75	0	if (distance < min_distance) {
76	0	min_distance = distance;
77	0	min = std::addressof(candidate);
78	0	}
79	0	}
80
81	0	if (min == nullptr) {
82	0	throw std::runtime_error("Failed to find next CoordinateChain");
83	0	}
84
85	0	return min;
86	0	}
87
88		template<typename T>
89		bool
90		hasMultipleUniqueCoordinates(const T& vec)
91	0	{
92	0	for (std::size_t i = 1; i < vec.size(); i++) {
93	0	if (vec[i] != vec[0]) {
94	0	return true;
95	0	}
96	0	}
97
98	0	return false;
99	0	}
100
101		/**
102		* @brief Identify counter-clockwise rings formed by the supplied coordinate vectors and the boundary of this box.
103		*
104		* @param box Box to be included in rings
105		* @param coord_lists A list of coordinate vectors, with each vector representing a traversal of `box` or a closed ring.
106		* @param visitor Function be applied to each ring identified. Because `coord_lists` may include both clockwise and
107		* counter-clockwise closed rings, `visitor` will be provided with the orientation of each ring as an
108		* argument.
109		*/
110		template<typename F>
111		void
112		visitRings(const Envelope& box, const std::vector<const std::vector<CoordinateXY>*>& coord_lists, F&& visitor)
113	0	{
114	0	std::vector<CoordinateChain> chains;
115	0	chains.reserve(coord_lists.size() + 4);
116
117	0	for (const auto& coords : coord_lists) {
118	0	if (!hasMultipleUniqueCoordinates(*coords)) {
119	0	continue;
120	0	}
121
122	0	if (coords->front() == coords->back() && hasMultipleUniqueCoordinates(*coords)) {
123		// Closed ring. Check orientation.
124
125		// TODO: Remove copy
126	0	CoordinateSequence seq(0, false, false);
127	0	seq.setPoints(*coords);
128	0	bool is_ccw = algorithm::Orientation::isCCW(&seq);
129	0	visitor(*coords, is_ccw);
130	0	} else {
131	0	double start = PerimeterDistance::getPerimeterDistance(box, coords->front());
132	0	double stop = PerimeterDistance::getPerimeterDistance(box, coords->back());
133
134	0	chains.emplace_back(start, stop, coords);
135	0	}
136	0	}
137
138	0	double height{ box.getHeight() };
139	0	double width{ box.getWidth() };
140	0	double perimeter{ box.getPerimeter() };
141
142		// create coordinate lists for corners
143	0	std::vector<CoordinateXY> bottom_left = { CoordinateXY(box.getMinX(), box.getMinY()) };
144	0	std::vector<CoordinateXY> top_left = { CoordinateXY(box.getMinX(), box.getMaxY()) };
145	0	std::vector<CoordinateXY> top_right = { CoordinateXY(box.getMaxX(), box.getMaxY()) };
146	0	std::vector<CoordinateXY> bottom_right = { CoordinateXY(box.getMaxX(), box.getMinY()) };
147
148		// Add chains for corners
149	0	chains.emplace_back(0.0, 0.0, &bottom_left);
150	0	chains.emplace_back(height, height, &top_left);
151	0	chains.emplace_back(height + width, height + width, &top_right);
152	0	chains.emplace_back(2 * height + width, 2 * height + width, &bottom_right);
153
154	0	std::vector<CoordinateXY> coords;
155	0	for (auto& chain_ref : chains) {
156	0	if (chain_ref.visited \|\| chain_ref.coordinates->size() == 1) {
157	0	continue;
158	0	}
159
160	0	coords.clear();
161
162	0	CoordinateChain* chain = std::addressof(chain_ref);
163	0	CoordinateChain* first_chain = chain;
164	0	do {
165	0	chain->visited = true;
166	0	coords.insert(coords.end(), chain->coordinates->cbegin(), chain->coordinates->cend());
167	0	chain = getNextChain(chains, chain, first_chain, perimeter);
168	0	} while (chain != first_chain);
169
170	0	coords.push_back(coords[0]);
171
172	0	if (hasMultipleUniqueCoordinates(coords)) {
173	0	visitor(coords, true);
174	0	}
175	0	}
176	0	} Unexecuted instantiation: TraversalAreas.cpp:void geos::operation::grid::visitRings<geos::operation::grid::TraversalAreas::getLeftHandArea(geos::geom::Envelope const&, std::__1::vector<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const, std::__1::allocator<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const> > const&)::$_0>(geos::geom::Envelope const&, std::__1::vector<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const, std::__1::allocator<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const> > const&, geos::operation::grid::TraversalAreas::getLeftHandArea(geos::geom::Envelope const&, std::__1::vector<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const, std::__1::allocator<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const> > const&)::$_0&&) Unexecuted instantiation: TraversalAreas.cpp:void geos::operation::grid::visitRings<geos::operation::grid::TraversalAreas::getLeftHandRings(geos::geom::GeometryFactory const&, geos::geom::Envelope const&, std::__1::vector<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const, std::__1::allocator<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const> > const&)::$_0>(geos::geom::Envelope const&, std::__1::vector<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const, std::__1::allocator<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const> > const&, geos::operation::grid::TraversalAreas::getLeftHandRings(geos::geom::GeometryFactory const&, geos::geom::Envelope const&, std::__1::vector<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const, std::__1::allocator<std::__1::vector<geos::geom::CoordinateXY, std::__1::allocator<geos::geom::CoordinateXY> > const> > const&)::$_0&&)
177
178		double
179		TraversalAreas::getLeftHandArea(const Envelope& box, const std::vector<const std::vector<CoordinateXY>*>& coord_lists)
180	0	{
181	0	double ccw_sum = 0;
182	0	double cw_sum = 0;
183	0	bool found_a_ring = false;
184
185	0	visitRings(box, coord_lists, [&cw_sum, &ccw_sum, &found_a_ring](const std::vector<CoordinateXY>& coords, bool is_ccw) {
186	0	found_a_ring = true;
187
188	0	if (is_ccw) {
189	0	ccw_sum += algorithm::Area::ofRing(coords);
190	0	} else {
191	0	cw_sum += algorithm::Area::ofRing(coords);
192	0	}
193	0	});
194
195	0	if (!found_a_ring) {
196	0	throw std::runtime_error("Cannot determine coverage fraction (it is either 0 or 100%)");
197	0	}
198
199		// If this box has only clockwise rings (holes) then the area
200		// of those holes should be subtracted from the complete box area.
201	0	if (ccw_sum < cw_sum) {
202	0	ccw_sum += box.getArea();
203	0	}
204
205	0	return ccw_sum - cw_sum;
206	0	}
207
208		std::unique_ptr<Geometry>
209		TraversalAreas::getLeftHandRings(const GeometryFactory& gfact, const Envelope& box, const std::vector<const std::vector<CoordinateXY>*>& coord_lists)
210	0	{
211	0	using geom::LinearRing;
212
213	0	std::vector<std::unique_ptr<LinearRing>> shells;
214	0	std::vector<std::unique_ptr<LinearRing>> holes;
215
216	0	bool found_a_ring = false;
217
218	0	visitRings(box, coord_lists, [&gfact, &shells, &holes, &found_a_ring](const std::vector<CoordinateXY>& coords, bool is_ccw) {
219	0	found_a_ring = true;
220
221		// finding a collapsed ring is sufficient to determine whether the cell interior is inside our outside,
222		// but we don't want to actually construct the ring.
223	0	if (algorithm::Area::ofRing(coords) == 0) {
224	0	return;
225	0	}
226
227	0	CoordinateSequence seq(0, false, false);
228	0	seq.setPoints(coords);
229	0	auto ring = gfact.createLinearRing(std::move(seq));
230
231	0	if (is_ccw) {
232	0	shells.push_back(std::move(ring));
233	0	} else {
234	0	holes.push_back(std::move(ring));
235	0	}
236	0	});
237
238	0	if (!found_a_ring) {
239	0	throw std::runtime_error("Cannot determine coverage fraction (it is either 0 or 100%)");
240	0	}
241
242	0	if (shells.empty() && holes.empty()) {
243	0	return gfact.createPolygon();
244	0	}
245
246	0	if (shells.empty() && !holes.empty()) {
247	0	CoordinateSequence seq(5, false, false);
248	0	seq.setAt(CoordinateXY{box.getMinX(), box.getMinY()}, 0);
249	0	seq.setAt(CoordinateXY{box.getMaxX(), box.getMinY()}, 1);
250	0	seq.setAt(CoordinateXY{box.getMaxX(), box.getMaxY()}, 2);
251	0	seq.setAt(CoordinateXY{box.getMinX(), box.getMaxY()}, 3);
252	0	seq.setAt(CoordinateXY{box.getMinX(), box.getMinY()}, 4);
253
254	0	shells.push_back(gfact.createLinearRing(std::move(seq)));
255	0	}
256
257	0	if (holes.empty()) {
258	0	std::vector<std::unique_ptr<Geometry>> polygons;
259	0	for (auto& shell : shells) {
260	0	polygons.push_back(gfact.createPolygon(std::move(shell)));
261	0	}
262
263	0	if (polygons.size() == 1) {
264	0	return std::move(polygons.front());
265	0	} else {
266	0	return gfact.createMultiPolygon(std::move(polygons));
267	0	}
268	0	}
269
270	0	if (shells.size() == 1) {
271	0	return gfact.createPolygon(std::move(shells.front()), std::move(holes));
272	0	}
273
274	0	polygonize::Polygonizer polygonizer(true);
275
276	0	for (const auto& shell : shells) {
277	0	polygonizer.add(static_cast<const Geometry*>(shell.get()));
278	0	}
279	0	for (const auto& hole : holes) {
280	0	polygonizer.add(static_cast<const Geometry*>(hole.get()));
281	0	}
282
283	0	auto polygonized = polygonizer.getPolygons();
284	0	return gfact.createMultiPolygon(std::move(polygonized));
285	0	}
286
287		}