-
Notifications
You must be signed in to change notification settings - Fork 310
/
OverlapUnion.cs
317 lines (278 loc) · 12.2 KB
/
OverlapUnion.cs
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
using System;
using System.Collections.Generic;
using NetTopologySuite.Geometries;
using NetTopologySuite.Geometries.Utilities;
namespace NetTopologySuite.Operation.Union
{
/// <summary>
/// Unions MultiPolygons efficiently by
/// using full topological union only for polygons which may overlap,
/// and combining with the remaining polygons.
/// Polygons which may overlap are those which intersect the common extent of the inputs.
/// Polygons wholly outside this extent must be disjoint to the computed union.
/// They can thus be simply combined with the union result,
/// which is much more performant.
/// (There is one caveat to this, which is discussed below).
/// <para/>
/// This situation is likely to occur during cascaded polygon union,
/// since the partitioning of polygons is done heuristically
/// and thus may group disjoint polygons which can lie far apart.
/// It may also occur in real world data which contains many disjoint polygons
/// (e.g. polygons representing parcels on different street blocks).
/// </summary>
/// <remarks>
/// <h2>Algorithm</h2>
/// The overlap region is determined as the common envelope of intersection.
/// The input polygons are partitioned into two sets:
/// <list type="bullet">
/// <item><term>Overlapping</term><description>Polygons which intersect the overlap region, and thus potentially overlap each other</description></item>
/// <item><term>Disjoint</term><description>Polygons which are disjoint from (lie wholly outside) the overlap region</description></item>
/// </list>
/// The Overlapping set is fully unioned, and then combined with the Disjoint set.
/// Performing a simple combine works because
/// the disjoint polygons do not interact with each
/// other(since the inputs are valid MultiPolygons).
/// They also do not interact with the Overlapping polygons,
/// since they are outside their envelope.
/// <h2>Discussion</h2>
/// In general the Overlapping set of polygons will
/// extend beyond the overlap envelope. This means that the union result
/// will extend beyond the overlap region.
/// There is a small chance that the topological
/// union of the overlap region will shift the result linework enough
/// that the result geometry intersects one of the Disjoint geometries.
/// This situation is detected and if it occurs
/// is remedied by falling back to performing a full union of the original inputs.
/// Detection is done by a fairly efficient comparison of edge segments which
/// extend beyond the overlap region. If any segments have changed
/// then there is a risk of introduced intersections, and full union is performed.
/// <para/>
/// This situation has not been observed in JTS using floating precision,
/// but it could happen due to snapping. It has been observed
/// in other APIs(e.g.GEOS) due to more aggressive snapping.
/// It is more likely to happen if a Snap - Rounding overlay is used.
/// <para/>
/// <b>NOTE: Test has shown that using this heuristic impairs performance.</b>
/// </remarks>
/// <author>Martin Davis</author>
[Obsolete("Due to impairing performance")]
public class OverlapUnion
{
/// <summary>
/// Union a pair of geometries,
/// using the more performant overlap union algorithm if possible.
/// </summary>
/// <param name="g0">A geometry to union</param>
/// <param name="g1">A geometry to union</param>
/// <returns>The union of the inputs</returns>
public static Geometry Union(Geometry g0, Geometry g1)
{
var union = new OverlapUnion(g0, g1);
return union.Union();
}
/// <summary>
/// Union a pair of geometries,
/// using the more performant overlap union algorithm if possible.
/// </summary>
/// <param name="g0">A geometry to union</param>
/// <param name="g1">A geometry to union</param>
/// <param name="unionFun">Function to union two geometries</param>
/// <returns>The union of the inputs</returns>
public static Geometry Union(Geometry g0, Geometry g1, UnionStrategy unionFun)
{
var union = new OverlapUnion(g0, g1, unionFun);
return union.Union();
}
private readonly GeometryFactory _geomFactory;
private readonly Geometry _g0;
private readonly Geometry _g1;
private readonly UnionStrategy _unionFun;
/// <summary>
/// Creates a new instance for unioning the given geometries.
/// </summary>
/// <param name="g0">A geometry to union</param>
/// <param name="g1">A geometry to union</param>
public OverlapUnion(Geometry g0, Geometry g1) : this(g0, g1, CascadedPolygonUnion.ClassicUnion)
{ }
/// <summary>
/// Creates a new instance for unioning the given geometries.
/// </summary>
/// <param name="g0">A geometry to union</param>
/// <param name="g1">A geometry to union</param>
/// <param name="unionFun">Function to union two geometries</param>
public OverlapUnion(Geometry g0, Geometry g1, UnionStrategy unionFun)
{
if (g0 == null)
{
throw new ArgumentNullException(nameof(g0));
}
_g0 = g0;
_g1 = g1;
_geomFactory = g0.Factory;
_unionFun = unionFun;
}
/// <summary>
/// Union a pair of geometries,
/// using the more performant overlap union algorithm if possible.
/// </summary>
/// <returns>The union of the inputs</returns>
public Geometry Union()
{
var overlapEnv = OverlapEnvelope(_g0, _g1);
/*
* If no overlap, can just combine the geometries
*/
if (overlapEnv.IsNull)
{
var g0Copy = _g0.Copy();
var g1Copy = _g1.Copy();
return GeometryCombiner.Combine(g0Copy, g1Copy);
}
var disjointPolys = new List<Geometry>();
var g0Overlap = ExtractByEnvelope(overlapEnv, _g0, disjointPolys);
var g1Overlap = ExtractByEnvelope(overlapEnv, _g1, disjointPolys);
//Console.WriteLine($"# geoms in common: {intersectingPolys.Count}");
var unionGeom = UnionFull(g0Overlap, g1Overlap);
Geometry result;
IsUnionOptimized = IsBorderSegmentsSame(unionGeom, overlapEnv);
if (!IsUnionOptimized)
{
// overlap union changed border segments... need to do full union
//System.out.println("OverlapUnion: Falling back to full union");
result = UnionFull(_g0, _g1);
}
else
{
//System.out.println("OverlapUnion: fast path");
result = Combine(unionGeom, disjointPolys);
}
return result;
}
/// <summary>
/// Gets a value indicating whether the optimized
/// or full union was performed.
/// </summary>
/// <remarks>Used for unit testing.</remarks>>
/// <returns><c>true</c> if the optimized union was performed</returns>
internal bool IsUnionOptimized { get; private set; }
private static Envelope OverlapEnvelope(Geometry g0, Geometry g1)
{
var g0Env = g0.EnvelopeInternal;
var g1Env = g1.EnvelopeInternal;
var overlapEnv = g0Env.Intersection(g1Env);
return overlapEnv;
}
private static Geometry Combine(Geometry unionGeom, List<Geometry> disjointPolys)
{
if (disjointPolys.Count <= 0)
return unionGeom;
disjointPolys.Add(unionGeom);
var result = GeometryCombiner.Combine(disjointPolys);
return result;
}
private Geometry ExtractByEnvelope(Envelope env, Geometry geom,
IList<Geometry> disjointGeoms)
{
var intersectingGeoms = new List<Geometry>();
for (int i = 0; i < geom.NumGeometries; i++)
{
var elem = geom.GetGeometryN(i);
if (elem.EnvelopeInternal.Intersects(env))
{
intersectingGeoms.Add(elem);
}
else
{
var copy = elem.Copy();
disjointGeoms.Add(copy);
}
}
return _geomFactory.BuildGeometry(intersectingGeoms);
}
private Geometry UnionFull(Geometry geom0, Geometry geom1)
{
// if both are empty collections, just return a copy of one of them
if (geom0.NumGeometries == 0
&& geom1.NumGeometries == 0)
return geom0.Copy();
var union = _unionFun.Union(geom0, geom1);
return union;
}
private bool IsBorderSegmentsSame(Geometry result, Envelope env)
{
var segsBefore = ExtractBorderSegments(_g0, _g1, env);
var segsAfter = new List<LineSegment>();
ExtractBorderSegments(result, env, segsAfter);
//System.out.println("# seg before: " + segsBefore.size() + " - # seg after: " + segsAfter.size());
return IsEqual(segsBefore, segsAfter);
}
private static bool IsEqual(ICollection<LineSegment> segs0, ICollection<LineSegment> segs1)
{
if (segs0.Count != segs1.Count)
return false;
var segIndex = new HashSet<LineSegment>(segs0);
foreach (var seg in segs1)
{
if (!segIndex.Contains(seg))
{
//System.out.println("Found changed border seg: " + seg);
return false;
}
}
return true;
}
private IList<LineSegment> ExtractBorderSegments(Geometry geom0, Geometry geom1, Envelope env)
{
var segs = new List<LineSegment>();
ExtractBorderSegments(geom0, env, segs);
if (geom1 != null)
ExtractBorderSegments(geom1, env, segs);
return segs;
}
private static void ExtractBorderSegments(Geometry geom, Envelope env, ICollection<LineSegment> segs)
{
geom.Apply(new BorderSegmentCoordinateFilter(env, segs));
}
private class BorderSegmentCoordinateFilter : ICoordinateSequenceFilter
{
private readonly ICollection<LineSegment> _segments;
private readonly Envelope _envelope;
public BorderSegmentCoordinateFilter(Envelope env, ICollection<LineSegment> segments)
{
_envelope = env;
_segments = segments;
}
public void Filter(CoordinateSequence seq, int i)
{
if (i <= 0) return;
// extract LineSegment
var p0 = seq.GetCoordinate(i - 1);
var p1 = seq.GetCoordinate(i);
bool isBorder = Intersects(_envelope, p0, p1) && !ContainsProperly(_envelope, p0, p1);
if (isBorder)
{
var seg = new LineSegment(p0, p1);
_segments.Add(seg);
}
}
public bool Done { get; } = false;
public bool GeometryChanged { get; } = false;
private static bool Intersects(Envelope env, Coordinate p0, Coordinate p1)
{
return env.Intersects(p0) || env.Intersects(p1);
}
private static bool ContainsProperly(Envelope env, Coordinate p0, Coordinate p1)
{
return ContainsProperly(env, p0) && ContainsProperly(env, p1);
}
private static bool ContainsProperly(Envelope env, Coordinate p)
{
if (env.IsNull) return false;
return p.X > env.MinX &&
p.X < env.MaxX &&
p.Y > env.MinY &&
p.Y < env.MaxY;
}
}
}
}