/src/geos/src/operation/overlayng/OverlayUtil.cpp

Source (jump to first uncovered line)
/**********************************************************************
 *
 * GEOS - Geometry Engine Open Source
 * http://geos.osgeo.org
 *
 * Copyright (C) 2020 Paul Ramsey <pramsey@cleverelephant.ca>
 *
 * 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 <geos/operation/overlayng/OverlayUtil.h>

#include <geos/operation/overlayng/InputGeometry.h>
#include <geos/operation/overlayng/OverlayGraph.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/PrecisionModel.h>
#include <geos/operation/overlayng/OverlayNG.h>
#include <geos/operation/overlayng/RobustClipEnvelopeComputer.h>
#include <geos/util/Assert.h>



namespace geos {      // geos
namespace operation { // geos.operation
namespace overlayng { // geos.operation.overlayng

using namespace geos::geom;

/*public static*/
bool
OverlayUtil::isFloating(const PrecisionModel* pm)
{
    if (pm == nullptr) return true;
    return pm->isFloating();
}

/*private static*/
double
OverlayUtil::safeExpandDistance(const Envelope* env, const PrecisionModel* pm)
{
    double envExpandDist;
    if (isFloating(pm)) {
        // if PM is FLOAT then there is no scale factor, so add 10%
        double minSize = std::min(env->getHeight(), env->getWidth());
        // heuristic to ensure zero-width envelopes don't cause total clipping
        if (minSize <= 0.0) {
            minSize = std::max(env->getHeight(), env->getWidth());
        }
        envExpandDist = SAFE_ENV_BUFFER_FACTOR * minSize;
    }
    else {
        // if PM is fixed, add a small multiple of the grid size
        double gridSize = 1.0 / pm->getScale();
        envExpandDist = SAFE_ENV_GRID_FACTOR * gridSize;
    }
    return envExpandDist;
}

/*private static*/
bool
OverlayUtil::safeEnv(const Envelope* env, const PrecisionModel* pm, Envelope& rsltEnvelope)
{
    double envExpandDist = safeExpandDistance(env, pm);
    rsltEnvelope = *env;
    rsltEnvelope.expandBy(envExpandDist);
    return true;
}

/*private static*/
bool
OverlayUtil::resultEnvelope(int opCode, const InputGeometry* inputGeom, const PrecisionModel* pm, Envelope& rsltEnvelope)
{
    switch (opCode) {
        case OverlayNG::INTERSECTION: {
            // use safe envelopes for intersection to ensure they contain rounded coordinates
            Envelope envA, envB;
            safeEnv(inputGeom->getEnvelope(0), pm, envA);
            safeEnv(inputGeom->getEnvelope(1), pm, envB);
            envA.intersection(envB, rsltEnvelope);
            return true;
        }
        case OverlayNG::DIFFERENCE: {
            safeEnv(inputGeom->getEnvelope(0), pm, rsltEnvelope);
            return true;
        }
    }
    // return false for UNION and SYMDIFFERENCE to indicate no clipping
    return false;
}

/*public static*/
bool
OverlayUtil::clippingEnvelope(int opCode, const InputGeometry* inputGeom, const PrecisionModel* pm, Envelope& rsltEnvelope)
{
    bool resultEnv = resultEnvelope(opCode, inputGeom, pm, rsltEnvelope);
    if (!resultEnv)
      return false;

    Envelope clipEnv = RobustClipEnvelopeComputer::getEnvelope(
        inputGeom->getGeometry(0),
        inputGeom->getGeometry(1),
        &rsltEnvelope);

    return safeEnv(&clipEnv, pm, rsltEnvelope);
}


/*public static*/
bool
OverlayUtil::isEmptyResult(int opCode, const Geometry* a, const Geometry* b, const PrecisionModel* pm)
{
    switch (opCode) {
    case OverlayNG::INTERSECTION: {
        if (isEnvDisjoint(a, b, pm))
            return true;
        break;
        }
    case OverlayNG::DIFFERENCE: {
        if (isEmpty(a))
            return true;
        break;
        }
    case OverlayNG::UNION:
    case OverlayNG::SYMDIFFERENCE: {
        if (isEmpty(a) && isEmpty(b))
            return true;
        break;
        }
    }
    return false;
}

/*private*/
bool
OverlayUtil::isEmpty(const Geometry* geom)
{
    return geom == nullptr || geom->isEmpty();
}

/*public static*/
bool
OverlayUtil::isEnvDisjoint(const Geometry* a, const Geometry* b, const PrecisionModel* pm)
{
    if (isEmpty(a) || isEmpty(b)) {
        return true;
    }
    if (isFloating(pm)) {
        return a->getEnvelopeInternal()->disjoint(b->getEnvelopeInternal());
    }
    return isDisjoint(a->getEnvelopeInternal(), b->getEnvelopeInternal(), pm);
}

/*private static*/
bool
OverlayUtil::isDisjoint(const Envelope* envA, const Envelope* envB, const PrecisionModel* pm)
{
    if (pm->makePrecise(envB->getMinX()) > pm->makePrecise(envA->getMaxX()))
        return true;
    if (pm->makePrecise(envB->getMaxX()) < pm->makePrecise(envA->getMinX()))
        return true;
    if (pm->makePrecise(envB->getMinY()) > pm->makePrecise(envA->getMaxY()))
        return true;
    if (pm->makePrecise(envB->getMaxY()) < pm->makePrecise(envA->getMinY()))
        return true;
    return false;
}

/*public static*/
std::unique_ptr<Geometry>
OverlayUtil::createEmptyResult(int dim, const GeometryFactory* geomFact)
{
    std::unique_ptr<Geometry> result(nullptr);
    switch (dim) {
    case 0:
        result = geomFact->createPoint();
        break;
    case 1:
        result = geomFact->createLineString();
        break;
    case 2:
        result = geomFact->createPolygon();
        break;
    case -1:
        result = geomFact->createGeometryCollection();
        break;
    default:
        util::Assert::shouldNeverReachHere("Unable to determine overlay result geometry dimension");
    }
    return result;
}

/*public static*/
int
OverlayUtil::resultDimension(int opCode, int dim0, int dim1)
{
    int resultDimension = -1;
    switch (opCode) {
    case OverlayNG::INTERSECTION:
        resultDimension = std::min(dim0, dim1);
        break;
    case OverlayNG::UNION:
        resultDimension = std::max(dim0, dim1);
        break;
    case OverlayNG::DIFFERENCE:
        resultDimension = dim0;
        break;
    case OverlayNG::SYMDIFFERENCE:
        /**
        * This result is chosen because
        * <pre>
        * SymDiff = Union( Diff(A, B), Diff(B, A) )
        * </pre>
        * and Union has the dimension of the highest-dimension argument.
        */
        resultDimension = std::max(dim0, dim1);
        break;
    }
    return resultDimension;
}

/* public static */
bool
OverlayUtil::isResultAreaConsistent(
    const Geometry* geom0, const Geometry* geom1,
    int opCode, const Geometry* result)
{
    if (geom0 == nullptr || geom1 == nullptr)
        return true;

    double areaResult = result->getArea();
    double areaA = geom0->getArea();
    double areaB = geom1->getArea();
    bool isConsistent = true;

    switch (opCode) {
    case OverlayNG::INTERSECTION:
        isConsistent = isLess(areaResult, areaA, AREA_HEURISTIC_TOLERANCE)
                    && isLess(areaResult, areaB, AREA_HEURISTIC_TOLERANCE);
        break;
    case OverlayNG::DIFFERENCE:
        isConsistent = isLess(areaResult, areaA, AREA_HEURISTIC_TOLERANCE)
                    && isGreater(areaResult, areaA - areaB, AREA_HEURISTIC_TOLERANCE);
        break;
    case OverlayNG::SYMDIFFERENCE:
        isConsistent = isLess(areaResult, areaA + areaB, AREA_HEURISTIC_TOLERANCE);
        break;
    case OverlayNG::UNION:
        isConsistent = isLess(areaA, areaResult, AREA_HEURISTIC_TOLERANCE)
                    && isLess(areaB, areaResult, AREA_HEURISTIC_TOLERANCE)
                    && isGreater(areaResult, areaA - areaB, AREA_HEURISTIC_TOLERANCE);
        break;
    }
    return isConsistent;
}


/*public static*/
std::unique_ptr<Geometry>
OverlayUtil::createResultGeometry(
    std::vector<std::unique_ptr<Polygon>>& resultPolyList,
    std::vector<std::unique_ptr<LineString>>& resultLineList,
    std::vector<std::unique_ptr<Point>>& resultPointList,
    const GeometryFactory* geometryFactory)
{
    std::vector<std::unique_ptr<Geometry>> geomList;

    // TODO: for mixed dimension, return collection of Multigeom for each dimension (breaking change)

    // element geometries of the result are always in the order A,L,P
    if (resultPolyList.size() > 0)
        moveGeometry(resultPolyList, geomList);
    if (resultLineList.size() > 0)
        moveGeometry(resultLineList, geomList);
    if (resultPointList.size() > 0)
        moveGeometry(resultPointList, geomList);

    // build the most specific geometry possible
    // TODO: perhaps do this internally to give more control?
    return geometryFactory->buildGeometry(std::move(geomList));
}

/*public static*/
std::unique_ptr<Geometry>
OverlayUtil::toLines(OverlayGraph* graph, bool isOutputEdges, const GeometryFactory* geomFact)
{
    std::vector<std::unique_ptr<LineString>> lines;
    std::vector<OverlayEdge*>& edges = graph->getEdges();
    for (OverlayEdge* edge : edges) {
      bool includeEdge = isOutputEdges || edge->isInResultArea();
      if (! includeEdge) continue;

      std::unique_ptr<CoordinateSequence> pts = edge->getCoordinatesOriented();
      std::unique_ptr<LineString> line = geomFact->createLineString(std::move(pts));
      // line->setUserData(labelForResult(edge));
      lines.push_back(std::move(line));
    }
    return geomFact->buildGeometry(std::move(lines));
}

/*public static*/
bool
OverlayUtil::round(const Point* pt, const PrecisionModel* pm, Coordinate& rsltCoord)
{
    if (pt->isEmpty()) return false;
    const Coordinate* p = pt->getCoordinate();
    rsltCoord = *p;
    if (! isFloating(pm)) {
        pm->makePrecise(rsltCoord);
    }
    return true;
}



} // namespace geos.operation.overlayng
} // namespace geos.operation
} // namespace geos

Coverage Report

Created: 2022-08-24 06:40

Line	Count	Source (jump to first uncovered line)
1		/**********************************************************************
2		*
3		* GEOS - Geometry Engine Open Source
4		* http://geos.osgeo.org
5		*
6		* Copyright (C) 2020 Paul Ramsey <pramsey@cleverelephant.ca>
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 <geos/operation/overlayng/OverlayUtil.h>
16
17		#include <geos/operation/overlayng/InputGeometry.h>
18		#include <geos/operation/overlayng/OverlayGraph.h>
19		#include <geos/geom/Coordinate.h>
20		#include <geos/geom/Envelope.h>
21		#include <geos/geom/GeometryFactory.h>
22		#include <geos/geom/PrecisionModel.h>
23		#include <geos/operation/overlayng/OverlayNG.h>
24		#include <geos/operation/overlayng/RobustClipEnvelopeComputer.h>
25		#include <geos/util/Assert.h>
26
27
28
29		namespace geos { // geos
30		namespace operation { // geos.operation
31		namespace overlayng { // geos.operation.overlayng
32
33		using namespace geos::geom;
34
35		/public static/
36		bool
37		OverlayUtil::isFloating(const PrecisionModel* pm)
38	251k	{
39	251k	if (pm == nullptr) return true;
40	216k	return pm->isFloating();
41	251k	}
42
43		/private static/
44		double
45		OverlayUtil::safeExpandDistance(const Envelope* env, const PrecisionModel* pm)
46	24.4k	{
47	24.4k	double envExpandDist;
48	24.4k	if (isFloating(pm)) {
49		// if PM is FLOAT then there is no scale factor, so add 10%
50	23.1k	double minSize = std::min(env->getHeight(), env->getWidth());
51		// heuristic to ensure zero-width envelopes don't cause total clipping
52	23.1k	if (minSize <= 0.0) {
53	912	minSize = std::max(env->getHeight(), env->getWidth());
54	912	}
55	23.1k	envExpandDist = SAFE_ENV_BUFFER_FACTOR * minSize;
56	23.1k	}
57	1.32k	else {
58		// if PM is fixed, add a small multiple of the grid size
59	1.32k	double gridSize = 1.0 / pm->getScale();
60	1.32k	envExpandDist = SAFE_ENV_GRID_FACTOR * gridSize;
61	1.32k	}
62	24.4k	return envExpandDist;
63	24.4k	}
64
65		/private static/
66		bool
67		OverlayUtil::safeEnv(const Envelope* env, const PrecisionModel* pm, Envelope& rsltEnvelope)
68	24.4k	{
69	24.4k	double envExpandDist = safeExpandDistance(env, pm);
70	24.4k	rsltEnvelope = *env;
71	24.4k	rsltEnvelope.expandBy(envExpandDist);
72	24.4k	return true;
73	24.4k	}
74
75		/private static/
76		bool
77		OverlayUtil::resultEnvelope(int opCode, const InputGeometry* inputGeom, const PrecisionModel* pm, Envelope& rsltEnvelope)
78	68.3k	{
79	68.3k	switch (opCode) {
80	4.79k	case OverlayNG::INTERSECTION: {
81		// use safe envelopes for intersection to ensure they contain rounded coordinates
82	4.79k	Envelope envA, envB;
83	4.79k	safeEnv(inputGeom->getEnvelope(0), pm, envA);
84	4.79k	safeEnv(inputGeom->getEnvelope(1), pm, envB);
85	4.79k	envA.intersection(envB, rsltEnvelope);
86	4.79k	return true;
87	0	}
88	5.04k	case OverlayNG::DIFFERENCE: {
89	5.04k	safeEnv(inputGeom->getEnvelope(0), pm, rsltEnvelope);
90	5.04k	return true;
91	0	}
92	68.3k	}
93		// return false for UNION and SYMDIFFERENCE to indicate no clipping
94	58.4k	return false;
95	68.3k	}
96
97		/public static/
98		bool
99		OverlayUtil::clippingEnvelope(int opCode, const InputGeometry* inputGeom, const PrecisionModel* pm, Envelope& rsltEnvelope)
100	68.3k	{
101	68.3k	bool resultEnv = resultEnvelope(opCode, inputGeom, pm, rsltEnvelope);
102	68.3k	if (!resultEnv)
103	58.4k	return false;
104
105	9.83k	Envelope clipEnv = RobustClipEnvelopeComputer::getEnvelope(
106	9.83k	inputGeom->getGeometry(0),
107	9.83k	inputGeom->getGeometry(1),
108	9.83k	&rsltEnvelope);
109
110	9.83k	return safeEnv(&clipEnv, pm, rsltEnvelope);
111	68.3k	}
112
113
114		/public static/
115		bool
116		OverlayUtil::isEmptyResult(int opCode, const Geometry* a, const Geometry* b, const PrecisionModel* pm)
117	124k	{
118	124k	switch (opCode) {
119	6.40k	case OverlayNG::INTERSECTION: {
120	6.40k	if (isEnvDisjoint(a, b, pm))
121	24	return true;
122	6.38k	break;
123	6.40k	}
124	8.98k	case OverlayNG::DIFFERENCE: {
125	8.98k	if (isEmpty(a))
126	38	return true;
127	8.95k	break;
128	8.98k	}
129	109k	case OverlayNG::UNION:
130	109k	case OverlayNG::SYMDIFFERENCE: {
131	109k	if (isEmpty(a) && isEmpty(b))
132	47.6k	return true;
133	61.7k	break;
134	109k	}
135	124k	}
136	77.0k	return false;
137	124k	}
138
139		/private/
140		bool
141		OverlayUtil::isEmpty(const Geometry* geom)
142	178k	{
143	178k	return geom == nullptr \|\| geom->isEmpty();
144	178k	}
145
146		/public static/
147		bool
148		OverlayUtil::isEnvDisjoint(const Geometry* a, const Geometry* b, const PrecisionModel* pm)
149	6.40k	{
150	6.40k	if (isEmpty(a) \|\| isEmpty(b)) {
151	13	return true;
152	13	}
153	6.39k	if (isFloating(pm)) {
154	6.13k	return a->getEnvelopeInternal()->disjoint(b->getEnvelopeInternal());
155	6.13k	}
156	265	return isDisjoint(a->getEnvelopeInternal(), b->getEnvelopeInternal(), pm);
157	6.39k	}
158
159		/private static/
160		bool
161		OverlayUtil::isDisjoint(const Envelope* envA, const Envelope* envB, const PrecisionModel* pm)
162	265	{
163	265	if (pm->makePrecise(envB->getMinX()) > pm->makePrecise(envA->getMaxX()))
164	0	return true;
165	265	if (pm->makePrecise(envB->getMaxX()) < pm->makePrecise(envA->getMinX()))
166	2	return true;
167	263	if (pm->makePrecise(envB->getMinY()) > pm->makePrecise(envA->getMaxY()))
168	0	return true;
169	263	if (pm->makePrecise(envB->getMaxY()) < pm->makePrecise(envA->getMinY()))
170	9	return true;
171	254	return false;
172	263	}
173
174		/public static/
175		std::unique_ptr<Geometry>
176		OverlayUtil::createEmptyResult(int dim, const GeometryFactory* geomFact)
177	51.3k	{
178	51.3k	std::unique_ptr<Geometry> result(nullptr);
179	51.3k	switch (dim) {
180	897	case 0:
181	897	result = geomFact->createPoint();
182	897	break;
183	681	case 1:
184	681	result = geomFact->createLineString();
185	681	break;
186	49.7k	case 2:
187	49.7k	result = geomFact->createPolygon();
188	49.7k	break;
189	0	case -1:
190	0	result = geomFact->createGeometryCollection();
191	0	break;
192	0	default:
193	0	util::Assert::shouldNeverReachHere("Unable to determine overlay result geometry dimension");
194	51.3k	}
195	51.3k	return result;
196	51.3k	}
197
198		/public static/
199		int
200		OverlayUtil::resultDimension(int opCode, int dim0, int dim1)
201	59.9k	{
202	59.9k	int resultDimension = -1;
203	59.9k	switch (opCode) {
204	2.05k	case OverlayNG::INTERSECTION:
205	2.05k	resultDimension = std::min(dim0, dim1);
206	2.05k	break;
207	53.7k	case OverlayNG::UNION:
208	53.7k	resultDimension = std::max(dim0, dim1);
209	53.7k	break;
210	4.13k	case OverlayNG::DIFFERENCE:
211	4.13k	resultDimension = dim0;
212	4.13k	break;
213	0	case OverlayNG::SYMDIFFERENCE:
214		/**
215		* This result is chosen because
216		* <pre>
217		* SymDiff = Union( Diff(A, B), Diff(B, A) )
218		* </pre>
219		* and Union has the dimension of the highest-dimension argument.
220		*/
221	0	resultDimension = std::max(dim0, dim1);
222	0	break;
223	59.9k	}
224	59.9k	return resultDimension;
225	59.9k	}
226
227		/* public static */
228		bool
229		OverlayUtil::isResultAreaConsistent(
230		const Geometry* geom0, const Geometry* geom1,
231		int opCode, const Geometry* result)
232	11.9k	{
233	11.9k	if (geom0 == nullptr \|\| geom1 == nullptr)
234	6.97k	return true;
235
236	5.00k	double areaResult = result->getArea();
237	5.00k	double areaA = geom0->getArea();
238	5.00k	double areaB = geom1->getArea();
239	5.00k	bool isConsistent = true;
240
241	5.00k	switch (opCode) {
242	1.88k	case OverlayNG::INTERSECTION:
243	1.88k	isConsistent = isLess(areaResult, areaA, AREA_HEURISTIC_TOLERANCE)
244	1.88k	&& isLess(areaResult, areaB, AREA_HEURISTIC_TOLERANCE);
245	1.88k	break;
246	1.70k	case OverlayNG::DIFFERENCE:
247	1.70k	isConsistent = isLess(areaResult, areaA, AREA_HEURISTIC_TOLERANCE)
248	1.70k	&& isGreater(areaResult, areaA - areaB, AREA_HEURISTIC_TOLERANCE);
249	1.70k	break;
250	0	case OverlayNG::SYMDIFFERENCE:
251	0	isConsistent = isLess(areaResult, areaA + areaB, AREA_HEURISTIC_TOLERANCE);
252	0	break;
253	1.41k	case OverlayNG::UNION:
254	1.41k	isConsistent = isLess(areaA, areaResult, AREA_HEURISTIC_TOLERANCE)
255	1.41k	&& isLess(areaB, areaResult, AREA_HEURISTIC_TOLERANCE)
256	1.41k	&& isGreater(areaResult, areaA - areaB, AREA_HEURISTIC_TOLERANCE);
257	1.41k	break;
258	5.00k	}
259	5.00k	return isConsistent;
260	5.00k	}
261
262
263		/public static/
264		std::unique_ptr<Geometry>
265		OverlayUtil::createResultGeometry(
266		std::vector<std::unique_ptr<Polygon>>& resultPolyList,
267		std::vector<std::unique_ptr<LineString>>& resultLineList,
268		std::vector<std::unique_ptr<Point>>& resultPointList,
269		const GeometryFactory* geometryFactory)
270	45.2k	{
271	45.2k	std::vector<std::unique_ptr<Geometry>> geomList;
272
273		// TODO: for mixed dimension, return collection of Multigeom for each dimension (breaking change)
274
275		// element geometries of the result are always in the order A,L,P
276	45.2k	if (resultPolyList.size() > 0)
277	40.9k	moveGeometry(resultPolyList, geomList);
278	45.2k	if (resultLineList.size() > 0)
279	4.59k	moveGeometry(resultLineList, geomList);
280	45.2k	if (resultPointList.size() > 0)
281	293	moveGeometry(resultPointList, geomList);
282
283		// build the most specific geometry possible
284		// TODO: perhaps do this internally to give more control?
285	45.2k	return geometryFactory->buildGeometry(std::move(geomList));
286	45.2k	}
287
288		/public static/
289		std::unique_ptr<Geometry>
290		OverlayUtil::toLines(OverlayGraph* graph, bool isOutputEdges, const GeometryFactory* geomFact)
291	0	{
292	0	std::vector<std::unique_ptr<LineString>> lines;
293	0	std::vector<OverlayEdge*>& edges = graph->getEdges();
294	0	for (OverlayEdge* edge : edges) {
295	0	bool includeEdge = isOutputEdges \|\| edge->isInResultArea();
296	0	if (! includeEdge) continue;
297
298	0	std::unique_ptr<CoordinateSequence> pts = edge->getCoordinatesOriented();
299	0	std::unique_ptr<LineString> line = geomFact->createLineString(std::move(pts));
300		// line->setUserData(labelForResult(edge));
301	0	lines.push_back(std::move(line));
302	0	}
303	0	return geomFact->buildGeometry(std::move(lines));
304	0	}
305
306		/public static/
307		bool
308		OverlayUtil::round(const Point* pt, const PrecisionModel* pm, Coordinate& rsltCoord)
309	75.1k	{
310	75.1k	if (pt->isEmpty()) return false;
311	75.1k	const Coordinate* p = pt->getCoordinate();
312	75.1k	rsltCoord = *p;
313	75.1k	if (! isFloating(pm)) {
314	420	pm->makePrecise(rsltCoord);
315	420	}
316	75.1k	return true;
317	75.1k	}
318
319
320
321		} // namespace geos.operation.overlayng
322		} // namespace geos.operation
323		} // namespace geos