/src/geos/src/operation/grid/Cell.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 <numeric>

#include <geos/algorithm/Length.h>
#include <geos/operation/grid/Cell.h>
#include <geos/operation/grid/Crossing.h>
#include <geos/operation/grid/TraversalAreas.h>

#define DEBUG_CELL 0
#include <iomanip>

#include "geos/geom/GeometryFactory.h"

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

namespace geos::operation::grid {

// Report where we leave an Envelope e when moving from c1 to c2
static Crossing
crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2)
{
    // vertical line
    if (c1.x == c2.x) {
        if (c2.y >= e.getMaxY()) {
            return Crossing{ Side::TOP, c1.x, e.getMaxY() };
        } else if (c2.y <= e.getMinY()) {
            return Crossing{ Side::BOTTOM, c1.x, e.getMinY() };
        //} else if (c2.x == e.getMinX()) {
        //    return Crossing{ Side::LEFT, c2.x, c2.y };
        //} else if (c2.x == e.getMaxX()) {
        //    return Crossing{ Side::RIGHT, c2.x, c2.y };
        } else {
            throw std::runtime_error("Never get here.");
        }
    }

    // horizontal line
    if (c1.y == c2.y) {
        if (c2.x >= e.getMaxX()) {
            return Crossing{ Side::RIGHT, e.getMaxX(), c1.y };
        } else if (c2.x <= e.getMinX()) {
            return Crossing{ Side::LEFT, e.getMinX(), c1.y };
        //} else if (c2.y == e.getMinY()) {
        //    return Crossing{ Side::BOTTOM, c2.x, c2.y };
        //} else if (c2.y == e.getMaxY()) {
        //    return Crossing{ Side::TOP, c2.x, c2.y };
        } else {
            throw std::runtime_error("Never get here");
        }
    }

    const double m = std::abs((c2.y - c1.y) / (c2.x - c1.x));

    const bool up = c2.y > c1.y;
    const bool right = c2.x > c1.x;

    if (up) {

        if (right) {
            // 1st quadrant
            const double y2 = c1.y + m * (e.getMaxX() - c1.x);

            if (y2 < e.getMaxY()) {
                return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
            } else {
                double x2 = c1.x + (e.getMaxY() - c1.y) / m;
                return Crossing{ Side::TOP, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMaxY() };
            }
        } else {
            // 2nd quadrant
            const double y2 = c1.y + m * (c1.x - e.getMinX());

            if (y2 < e.getMaxY()) {
                return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
            } else {
                double x2 = c1.x - (e.getMaxY() - c1.y) / m;
                return Crossing{ Side::TOP, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMaxY() };
            }
        }
    } else {
        // For downward segments, we calculate constructed coordinates relative to c2, not c1. This is so the
        // same coordinate will be calculated regardless of the segment orientation. This is important for maintaining
        // valid polygon coverages (the same segment will be processed with opposite orientation along shared
        // boundaries)

        if (right) {
            // 4th quadrant
            const double y2 = c2.y + m * (c2.x - e.getMaxX());

            if (y2 > e.getMinY()) {
                return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
            } else {
                double x2 = c2.x - (e.getMinY() - c2.y) / m;
                return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() };
            }
        } else {
            // 3rd quadrant
            const double y2 = c2.y + m * (e.getMinX() - c2.x);

            if (y2 > e.getMinY()) {
                return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
            } else {
                double x2 = c2.x + (e.getMinY() - c2.y) / m;
                return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() };
            }
        }
    }
}


double
Cell::getHeight() const
{
    return m_box.getHeight();
}

double
Cell::getWidth() const
{
    return m_box.getWidth();
}

double
Cell::getArea() const
{
    return m_box.getArea();
}

Side
Cell::getSide(const CoordinateXY& c) const
{
    if (c.x == m_box.getMinX()) {
        return Side::LEFT;
    } else if (c.x == m_box.getMaxX()) {
        return Side::RIGHT;
    } else if (c.y == m_box.getMinY()) {
        return Side::BOTTOM;
    } else if (c.y == m_box.getMaxY()) {
        return Side::TOP;
    }

    return Side::NONE;
}

void
Cell::forceExit()
{
    if (getLastTraversal().isExited()) {
        return;
    }

    const CoordinateXY& last = getLastTraversal().getLastCoordinate();

    if (getLocation(last) == Location::BOUNDARY) {
        getLastTraversal().forceExit(getSide(last));
    }
}

Cell::Location
Cell::getLocation(const CoordinateXY& c) const
{
    if (m_box.containsProperly(c)) {
        return Cell::Location::INSIDE;
    }

    if (m_box.contains(c)) {
        return Cell::Location::BOUNDARY;
    }

    return Cell::Location::OUTSIDE;
}

Traversal&
Cell::traversal_in_progress()
{
    if (m_traversals.empty() || m_traversals.back().isExited() || m_traversals.back().isClosedRing()) {
        m_traversals.emplace_back();
    }

    return m_traversals[m_traversals.size() - 1];
}

Traversal&
Cell::getLastTraversal()
{
    return m_traversals.at(m_traversals.size() - 1);
}

bool
Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original, bool exitOnBoundary, const void* parentage)
{
    Traversal& t = traversal_in_progress();

    if (t.isEmpty()) {
#if DEBUG_CELL
        std::cout << std::setprecision(17);
        std::cout << "Entering " << m_box << " from " << getSide(c) << " at " << c << std::endl;
#endif

        t.enter(c, getSide(c), parentage);
        return true;
    }

    const auto loc = getLocation(c);
    const bool canTakeCoordinate = exitOnBoundary ? (loc == Location::INSIDE) : (loc != Location::OUTSIDE);

    if (canTakeCoordinate) {
#if DEBUG_CELL
        std::cout << "Still in " << m_box << " with " << c << std::endl;
#endif

        t.add(c);

        if (t.isClosedRing()) {
            t.forceExit(Side::NONE);
        }

        return true;
    } else if (loc == Location::BOUNDARY) {
        if (c.x == m_box.getMaxX()) {
            t.exit(c, Side::RIGHT);
        } else if (c.y == m_box.getMaxY()) {
            t.exit(c, Side::TOP);
        } else if (c.x == m_box.getMinX()) {
            t.exit(c, Side::LEFT);
        } else {
            t.exit(c, Side::BOTTOM);
        }
        return false;
    }

    // We need to calculate the coordinate of the cell exit point using only uninterpolated coordinates.
    // (The previous point in the traversal may be an interpolated coordinate.) If an interpolated coordinate
    // is used, it can cause an error in the relative position two traversals, inverting the fraction of
    // the cell that is considered covered. (See robustness regression test #7).
    Crossing x = prev_original ? crossing(m_box, *prev_original, c) : crossing(m_box, t.getLastCoordinate(), c);
    t.exit(x.getCoord(), x.getSide());

#if DEBUG_CELL
    std::cout << "Leaving " << m_box << " from " << x.getSide() << " at " << x.getCoord();
    std::cout << " on the way to " << c << std::endl;
#endif

    return false;
}

double
Cell::getTraversalLength() const
{
    return std::accumulate(m_traversals.begin(), m_traversals.end(), 0.0, [](double tot, const Traversal& t) {
        return tot + algorithm::Length::ofLine(t.getCoordinates());
    });
}

bool
Cell::isDetermined() const
{
    for (const auto& t : m_traversals) {
        if (t.isClosedRing()) {
            if (!t.isClosedRingWithArea()) {
                continue;
            }
        } else if (!t.isTraversed()) {
            continue;
        }

        if (t.hasMultipleUniqueCoordinates()) {
            return true;
        }
    }

    return false;
}

std::vector<const Traversal*>
Cell::getTraversals() const
{
    std::vector<const Traversal*> traversals;
    traversals.reserve(m_traversals.size());

    for (const auto& t : m_traversals) {
        if (t.isTraversed() || t.isClosedRingWithArea()) {
            traversals.push_back(&t);
        }
    }

    return traversals;
}

double
Cell::getCoveredFraction() const
{
    auto coord_lists = getTraversals();
    return TraversalAreas::getLeftHandArea(m_box, coord_lists) / getArea();
}

std::unique_ptr<Geometry>
Cell::getCoveredPolygons(const GeometryFactory& gfact) const
{
    auto coord_lists = getTraversals();
    return TraversalAreas::getLeftHandRings(gfact, m_box, coord_lists);
}

void Cell::getEdgePoints(Side s, std::vector<CoordinateXY> &edgePoints) const
{
    for (const Traversal& t : m_traversals) {
        const auto& coords = t.getCoordinates();

        for (const auto& c : coords) {
            if ((s == Side::LEFT && c.x == m_box.getMinX()) ||
                (s == Side::RIGHT && c.x == m_box.getMaxX()) ||
                (s == Side::BOTTOM && c.y == m_box.getMinY()) ||
                (s == Side::TOP && c.y == m_box.getMaxY())) {
                edgePoints.push_back(c);
            }
        }
    }
}

}

Coverage Report

Created: 2026-07-16 06:39

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 <numeric>
16
17		#include <geos/algorithm/Length.h>
18		#include <geos/operation/grid/Cell.h>
19		#include <geos/operation/grid/Crossing.h>
20		#include <geos/operation/grid/TraversalAreas.h>
21
22		#define DEBUG_CELL 0
23		#include <iomanip>
24
25		#include "geos/geom/GeometryFactory.h"
26
27		using geos::geom::CoordinateXY;
28		using geos::geom::Geometry;
29		using geos::geom::GeometryFactory;
30
31		namespace geos::operation::grid {
32
33		// Report where we leave an Envelope e when moving from c1 to c2
34		static Crossing
35		crossing(const geom::Envelope& e, const CoordinateXY& c1, const CoordinateXY& c2)
36	0	{
37		// vertical line
38	0	if (c1.x == c2.x) {
39	0	if (c2.y >= e.getMaxY()) {
40	0	return Crossing{ Side::TOP, c1.x, e.getMaxY() };
41	0	} else if (c2.y <= e.getMinY()) {
42	0	return Crossing{ Side::BOTTOM, c1.x, e.getMinY() };
43		//} else if (c2.x == e.getMinX()) {
44		// return Crossing{ Side::LEFT, c2.x, c2.y };
45		//} else if (c2.x == e.getMaxX()) {
46		// return Crossing{ Side::RIGHT, c2.x, c2.y };
47	0	} else {
48	0	throw std::runtime_error("Never get here.");
49	0	}
50	0	}
51
52		// horizontal line
53	0	if (c1.y == c2.y) {
54	0	if (c2.x >= e.getMaxX()) {
55	0	return Crossing{ Side::RIGHT, e.getMaxX(), c1.y };
56	0	} else if (c2.x <= e.getMinX()) {
57	0	return Crossing{ Side::LEFT, e.getMinX(), c1.y };
58		//} else if (c2.y == e.getMinY()) {
59		// return Crossing{ Side::BOTTOM, c2.x, c2.y };
60		//} else if (c2.y == e.getMaxY()) {
61		// return Crossing{ Side::TOP, c2.x, c2.y };
62	0	} else {
63	0	throw std::runtime_error("Never get here");
64	0	}
65	0	}
66
67	0	const double m = std::abs((c2.y - c1.y) / (c2.x - c1.x));
68
69	0	const bool up = c2.y > c1.y;
70	0	const bool right = c2.x > c1.x;
71
72	0	if (up) {
73
74	0	if (right) {
75		// 1st quadrant
76	0	const double y2 = c1.y + m * (e.getMaxX() - c1.x);
77
78	0	if (y2 < e.getMaxY()) {
79	0	return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
80	0	} else {
81	0	double x2 = c1.x + (e.getMaxY() - c1.y) / m;
82	0	return Crossing{ Side::TOP, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMaxY() };
83	0	}
84	0	} else {
85		// 2nd quadrant
86	0	const double y2 = c1.y + m * (c1.x - e.getMinX());
87
88	0	if (y2 < e.getMaxY()) {
89	0	return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
90	0	} else {
91	0	double x2 = c1.x - (e.getMaxY() - c1.y) / m;
92	0	return Crossing{ Side::TOP, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMaxY() };
93	0	}
94	0	}
95	0	} else {
96		// For downward segments, we calculate constructed coordinates relative to c2, not c1. This is so the
97		// same coordinate will be calculated regardless of the segment orientation. This is important for maintaining
98		// valid polygon coverages (the same segment will be processed with opposite orientation along shared
99		// boundaries)
100
101	0	if (right) {
102		// 4th quadrant
103	0	const double y2 = c2.y + m * (c2.x - e.getMaxX());
104
105	0	if (y2 > e.getMinY()) {
106	0	return Crossing{ Side::RIGHT, e.getMaxX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
107	0	} else {
108	0	double x2 = c2.x - (e.getMinY() - c2.y) / m;
109	0	return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() };
110	0	}
111	0	} else {
112		// 3rd quadrant
113	0	const double y2 = c2.y + m * (e.getMinX() - c2.x);
114
115	0	if (y2 > e.getMinY()) {
116	0	return Crossing{ Side::LEFT, e.getMinX(), std::clamp(y2, e.getMinY(), e.getMaxY()) };
117	0	} else {
118	0	double x2 = c2.x + (e.getMinY() - c2.y) / m;
119	0	return Crossing{ Side::BOTTOM, std::clamp(x2, e.getMinX(), e.getMaxX()), e.getMinY() };
120	0	}
121	0	}
122	0	}
123	0	}
124
125
126		double
127		Cell::getHeight() const
128	0	{
129	0	return m_box.getHeight();
130	0	}
131
132		double
133		Cell::getWidth() const
134	0	{
135	0	return m_box.getWidth();
136	0	}
137
138		double
139		Cell::getArea() const
140	0	{
141	0	return m_box.getArea();
142	0	}
143
144		Side
145		Cell::getSide(const CoordinateXY& c) const
146	0	{
147	0	if (c.x == m_box.getMinX()) {
148	0	return Side::LEFT;
149	0	} else if (c.x == m_box.getMaxX()) {
150	0	return Side::RIGHT;
151	0	} else if (c.y == m_box.getMinY()) {
152	0	return Side::BOTTOM;
153	0	} else if (c.y == m_box.getMaxY()) {
154	0	return Side::TOP;
155	0	}
156
157	0	return Side::NONE;
158	0	}
159
160		void
161		Cell::forceExit()
162	0	{
163	0	if (getLastTraversal().isExited()) {
164	0	return;
165	0	}
166
167	0	const CoordinateXY& last = getLastTraversal().getLastCoordinate();
168
169	0	if (getLocation(last) == Location::BOUNDARY) {
170	0	getLastTraversal().forceExit(getSide(last));
171	0	}
172	0	}
173
174		Cell::Location
175		Cell::getLocation(const CoordinateXY& c) const
176	0	{
177	0	if (m_box.containsProperly(c)) {
178	0	return Cell::Location::INSIDE;
179	0	}
180
181	0	if (m_box.contains(c)) {
182	0	return Cell::Location::BOUNDARY;
183	0	}
184
185	0	return Cell::Location::OUTSIDE;
186	0	}
187
188		Traversal&
189		Cell::traversal_in_progress()
190	0	{
191	0	if (m_traversals.empty() \|\| m_traversals.back().isExited() \|\| m_traversals.back().isClosedRing()) {
192	0	m_traversals.emplace_back();
193	0	}
194
195	0	return m_traversals[m_traversals.size() - 1];
196	0	}
197
198		Traversal&
199		Cell::getLastTraversal()
200	0	{
201	0	return m_traversals.at(m_traversals.size() - 1);
202	0	}
203
204		bool
205		Cell::take(const CoordinateXY& c, const CoordinateXY* prev_original, bool exitOnBoundary, const void* parentage)
206	0	{
207	0	Traversal& t = traversal_in_progress();
208
209	0	if (t.isEmpty()) {
210		#if DEBUG_CELL
211		std::cout << std::setprecision(17);
212		std::cout << "Entering " << m_box << " from " << getSide(c) << " at " << c << std::endl;
213		#endif
214
215	0	t.enter(c, getSide(c), parentage);
216	0	return true;
217	0	}
218
219	0	const auto loc = getLocation(c);
220	0	const bool canTakeCoordinate = exitOnBoundary ? (loc == Location::INSIDE) : (loc != Location::OUTSIDE);
221
222	0	if (canTakeCoordinate) {
223		#if DEBUG_CELL
224		std::cout << "Still in " << m_box << " with " << c << std::endl;
225		#endif
226
227	0	t.add(c);
228
229	0	if (t.isClosedRing()) {
230	0	t.forceExit(Side::NONE);
231	0	}
232
233	0	return true;
234	0	} else if (loc == Location::BOUNDARY) {
235	0	if (c.x == m_box.getMaxX()) {
236	0	t.exit(c, Side::RIGHT);
237	0	} else if (c.y == m_box.getMaxY()) {
238	0	t.exit(c, Side::TOP);
239	0	} else if (c.x == m_box.getMinX()) {
240	0	t.exit(c, Side::LEFT);
241	0	} else {
242	0	t.exit(c, Side::BOTTOM);
243	0	}
244	0	return false;
245	0	}
246
247		// We need to calculate the coordinate of the cell exit point using only uninterpolated coordinates.
248		// (The previous point in the traversal may be an interpolated coordinate.) If an interpolated coordinate
249		// is used, it can cause an error in the relative position two traversals, inverting the fraction of
250		// the cell that is considered covered. (See robustness regression test #7).
251	0	Crossing x = prev_original ? crossing(m_box, *prev_original, c) : crossing(m_box, t.getLastCoordinate(), c);
252	0	t.exit(x.getCoord(), x.getSide());
253
254		#if DEBUG_CELL
255		std::cout << "Leaving " << m_box << " from " << x.getSide() << " at " << x.getCoord();
256		std::cout << " on the way to " << c << std::endl;
257		#endif
258
259	0	return false;
260	0	}
261
262		double
263		Cell::getTraversalLength() const
264	0	{
265	0	return std::accumulate(m_traversals.begin(), m_traversals.end(), 0.0, [](double tot, const Traversal& t) {
266	0	return tot + algorithm::Length::ofLine(t.getCoordinates());
267	0	});
268	0	}
269
270		bool
271		Cell::isDetermined() const
272	0	{
273	0	for (const auto& t : m_traversals) {
274	0	if (t.isClosedRing()) {
275	0	if (!t.isClosedRingWithArea()) {
276	0	continue;
277	0	}
278	0	} else if (!t.isTraversed()) {
279	0	continue;
280	0	}
281
282	0	if (t.hasMultipleUniqueCoordinates()) {
283	0	return true;
284	0	}
285	0	}
286
287	0	return false;
288	0	}
289
290		std::vector<const Traversal*>
291		Cell::getTraversals() const
292	0	{
293	0	std::vector<const Traversal*> traversals;
294	0	traversals.reserve(m_traversals.size());
295
296	0	for (const auto& t : m_traversals) {
297	0	if (t.isTraversed() \|\| t.isClosedRingWithArea()) {
298	0	traversals.push_back(&t);
299	0	}
300	0	}
301
302	0	return traversals;
303	0	}
304
305		double
306		Cell::getCoveredFraction() const
307	0	{
308	0	auto coord_lists = getTraversals();
309	0	return TraversalAreas::getLeftHandArea(m_box, coord_lists) / getArea();
310	0	}
311
312		std::unique_ptr<Geometry>
313		Cell::getCoveredPolygons(const GeometryFactory& gfact) const
314	0	{
315	0	auto coord_lists = getTraversals();
316	0	return TraversalAreas::getLeftHandRings(gfact, m_box, coord_lists);
317	0	}
318
319		void Cell::getEdgePoints(Side s, std::vector<CoordinateXY> &edgePoints) const
320	0	{
321	0	for (const Traversal& t : m_traversals) {
322	0	const auto& coords = t.getCoordinates();
323
324	0	for (const auto& c : coords) {
325	0	if ((s == Side::LEFT && c.x == m_box.getMinX()) \|\|
326	0	(s == Side::RIGHT && c.x == m_box.getMaxX()) \|\|
327	0	(s == Side::BOTTOM && c.y == m_box.getMinY()) \|\|
328	0	(s == Side::TOP && c.y == m_box.getMaxY())) {
329	0	edgePoints.push_back(c);
330	0	}
331	0	}
332	0	}
333	0	}
334
335		}