1 /**********************************************************************
2 *
3 * GEOS - Geometry Engine Open Source
4 * http://geos.osgeo.org
5 *
6 * Copyright (C) 2009-2011 Sandro Santilli <strk@kbt.io>
7 * Copyright (C) 2005-2007 Refractions Research Inc.
8 * Copyright (C) 2001-2002 Vivid Solutions Inc.
9 *
10 * This is free software; you can redistribute and/or modify it under
11 * the terms of the GNU Lesser General Public Licence as published
12 * by the Free Software Foundation.
13 * See the COPYING file for more information.
14 *
15 **********************************************************************
16 *
17 * Last port: operation/buffer/BufferOp.java r378 (JTS-1.12)
18 *
19 **********************************************************************/
20
21 #include <algorithm>
22 #include <cmath>
23
24 #include <geos/constants.h>
25 #include <geos/profiler.h>
26 #include <geos/precision/GeometryPrecisionReducer.h>
27 #include <geos/operation/buffer/BufferOp.h>
28 #include <geos/operation/buffer/BufferBuilder.h>
29 #include <geos/geom/Geometry.h>
30 #include <geos/geom/GeometryFactory.h>
31 #include <geos/geom/PrecisionModel.h>
32
33 #include <geos/noding/ScaledNoder.h>
34
35 #include <geos/noding/snapround/MCIndexSnapRounder.h>
36 #include <geos/noding/snapround/MCIndexPointSnapper.h>
37 #include <geos/noding/snapround/SnapRoundingNoder.h>
38
39 //FIXME: for temporary use, see other FIXME in file
40 #include <geos/algorithm/LineIntersector.h>
41 #include <geos/noding/MCIndexNoder.h>
42 #include <geos/noding/IntersectionAdder.h>
43
44
45
46
47 #ifndef GEOS_DEBUG
48 #define GEOS_DEBUG 0
49 #endif
50
51 //#define PROFILE 1
52
53 using namespace geos::noding;
54 using namespace geos::geom;
55
56 namespace geos {
57 namespace operation { // geos.operation
58 namespace buffer { // geos.operation.buffer
59
60 #if PROFILE
61 static Profiler* profiler = Profiler::instance();
62 #endif
63
64 #if 0
65 double
66 OLDprecisionScaleFactor(const Geometry* g,
67 double distance, int maxPrecisionDigits)
68 {
69 const Envelope* env = g->getEnvelopeInternal();
70 double envSize = std::max(env->getHeight(), env->getWidth());
71 double expandByDistance = distance > 0.0 ? distance : 0.0;
72 double bufEnvSize = envSize + 2 * expandByDistance;
73 // the smallest power of 10 greater than the buffer envelope
74 int bufEnvLog10 = (int)(std::log(bufEnvSize) / std::log(10.0) + 1.0);
75 int minUnitLog10 = bufEnvLog10 - maxPrecisionDigits;
76 // scale factor is inverse of min Unit size, so flip sign of exponent
77 double scaleFactor = std::pow(10.0, -minUnitLog10);
78 return scaleFactor;
79 }
80 #endif
81
82 /*private*/
83 double
precisionScaleFactor(const Geometry * g,double distance,int maxPrecisionDigits)84 BufferOp::precisionScaleFactor(const Geometry* g,
85 double distance,
86 int maxPrecisionDigits)
87 {
88 const Envelope* env = g->getEnvelopeInternal();
89 double envMax = std::max(
90 std::max(fabs(env->getMaxX()), fabs(env->getMinX())),
91 std::max(fabs(env->getMaxY()), fabs(env->getMinY()))
92 );
93
94 double expandByDistance = distance > 0.0 ? distance : 0.0;
95 double bufEnvMax = envMax + 2 * expandByDistance;
96
97 // the smallest power of 10 greater than the buffer envelope
98 int bufEnvPrecisionDigits = (int)(std::log(bufEnvMax) / std::log(10.0) + 1.0);
99 int minUnitLog10 = maxPrecisionDigits - bufEnvPrecisionDigits;
100
101 double scaleFactor = std::pow(10.0, minUnitLog10);
102
103 return scaleFactor;
104 }
105
106 /*public static*/
107 Geometry*
bufferOp(const Geometry * g,double distance,int quadrantSegments,int nEndCapStyle)108 BufferOp::bufferOp(const Geometry* g, double distance,
109 int quadrantSegments,
110 int nEndCapStyle)
111 {
112 BufferOp bufOp(g);
113 bufOp.setQuadrantSegments(quadrantSegments);
114 bufOp.setEndCapStyle(nEndCapStyle);
115 return bufOp.getResultGeometry(distance);
116 }
117
118 /*public*/
119 Geometry*
getResultGeometry(double nDistance)120 BufferOp::getResultGeometry(double nDistance)
121 {
122 distance = nDistance;
123 computeGeometry();
124 return resultGeometry;
125 }
126
127 /*private*/
128 void
computeGeometry()129 BufferOp::computeGeometry()
130 {
131 #if GEOS_DEBUG
132 std::cerr << "BufferOp::computeGeometry: trying with original precision" << std::endl;
133 #endif
134
135 bufferOriginalPrecision();
136
137 if(resultGeometry != nullptr) {
138 return;
139 }
140
141 #if GEOS_DEBUG
142 std::cerr << "bufferOriginalPrecision failed (" << saveException.what() << "), trying with reduced precision"
143 << std::endl;
144 #endif
145
146 const PrecisionModel& argPM = *(argGeom->getFactory()->getPrecisionModel());
147 if(argPM.getType() == PrecisionModel::FIXED) {
148 bufferFixedPrecision(argPM);
149 }
150 else {
151 bufferReducedPrecision();
152 }
153 }
154
155 /*private*/
156 void
bufferReducedPrecision()157 BufferOp::bufferReducedPrecision()
158 {
159
160 // try and compute with decreasing precision,
161 // up to a min, to avoid gross results
162 // (not in JTS, see http://trac.osgeo.org/geos/ticket/605)
163 #define MIN_PRECISION_DIGITS 6
164 for(int precDigits = MAX_PRECISION_DIGITS; precDigits >= MIN_PRECISION_DIGITS; precDigits--) {
165 #if GEOS_DEBUG
166 std::cerr << "BufferOp::computeGeometry: trying with precDigits " << precDigits << std::endl;
167 #endif
168 try {
169 bufferReducedPrecision(precDigits);
170 }
171 catch(const util::TopologyException& ex) {
172 saveException = ex;
173 // don't propagate the exception - it will be detected by fact that resultGeometry is null
174 }
175
176 if(resultGeometry != nullptr) {
177 // debug
178 //if ( saveException ) std::cerr<<saveException->toString()<<std::endl;
179 return;
180 }
181 }
182 // tried everything - have to bail
183 throw saveException;
184 }
185
186 /*private*/
187 void
bufferOriginalPrecision()188 BufferOp::bufferOriginalPrecision()
189 {
190 BufferBuilder bufBuilder(bufParams);
191
192 //std::cerr<<"computing with original precision"<<std::endl;
193 try {
194 resultGeometry = bufBuilder.buffer(argGeom, distance);
195 }
196 catch(const util::TopologyException& ex) {
197 // don't propagate the exception - it will be detected by
198 // fact that resultGeometry is null
199 saveException = ex;
200
201 //std::cerr<<ex->toString()<<std::endl;
202 }
203 //std::cerr<<"done"<<std::endl;
204 }
205
206 void
bufferReducedPrecision(int precisionDigits)207 BufferOp::bufferReducedPrecision(int precisionDigits)
208 {
209 double sizeBasedScaleFactor = precisionScaleFactor(argGeom, distance, precisionDigits);
210
211 #if GEOS_DEBUG
212 std::cerr << "recomputing with precision scale factor = "
213 << sizeBasedScaleFactor
214 << std::endl;
215 #endif
216
217 assert(sizeBasedScaleFactor > 0);
218 PrecisionModel fixedPM(sizeBasedScaleFactor);
219 bufferFixedPrecision(fixedPM);
220 }
221
222 /*private*/
223 void
bufferFixedPrecision(const PrecisionModel & fixedPM)224 BufferOp::bufferFixedPrecision(const PrecisionModel& fixedPM)
225 {
226 PrecisionModel pm(1.0); // fixed as well
227
228 #define SNAP_WITH_NODER
229 #ifdef SNAP_WITH_NODER
230 // Reduce precision using SnapRoundingNoder
231 //
232 // This more closely aligns with JTS implementation,
233 // and avoids reducing the precision of the input
234 // geometry.
235 //
236 // TODO: Add a finer fallback sequence. Full
237 // precision, then SnappingNoder, then
238 // SnapRoundingNoder.
239
240 snapround::SnapRoundingNoder inoder(&pm);
241 ScaledNoder noder(inoder, fixedPM.getScale());
242 BufferBuilder bufBuilder(bufParams);
243 bufBuilder.setWorkingPrecisionModel(&fixedPM);
244 bufBuilder.setNoder(&noder);
245 resultGeometry = bufBuilder.buffer(argGeom, distance);
246
247 #else
248 algorithm::LineIntersector li(&fixedPM);
249 IntersectionAdder ia(li);
250 MCIndexNoder inoder(&ia);
251 ScaledNoder noder(inoder, fixedPM.getScale());
252 BufferBuilder bufBuilder(bufParams);
253 bufBuilder.setWorkingPrecisionModel(&fixedPM);
254 bufBuilder.setNoder(&noder);
255
256 // Snap by reducing the precision of the input geometry
257 //
258 // NOTE: this reduction is not in JTS and should supposedly
259 // not be needed because the PrecisionModel we pass
260 // to the BufferBuilder above (with setWorkingPrecisionModel)
261 // should be used to round coordinates emitted by the
262 // OffsetCurveBuilder, thus effectively producing a fully
263 // rounded input to the noder.
264 // Nonetheless the amount of scrambling done by rounding here
265 // is known to fix at least one case in which MCIndexNoder
266 // would fail: http://trac.osgeo.org/geos/ticket/605
267 //
268 // TODO: follow JTS in MCIndexSnapRounder usage
269 //
270 const Geometry* workGeom = argGeom;
271 const PrecisionModel& argPM = *(argGeom->getFactory()->getPrecisionModel());
272 std::unique_ptr<Geometry> fixedGeom;
273 if(argPM.getType() != PrecisionModel::FIXED || argPM.getScale() != fixedPM.getScale()) {
274 using precision::GeometryPrecisionReducer;
275 fixedGeom = GeometryPrecisionReducer::reduce(*argGeom, fixedPM);
276 workGeom = fixedGeom.get();
277 }
278
279 // this may throw an exception, if robustness errors are encountered
280 resultGeometry = bufBuilder.buffer(workGeom, distance);
281 #endif
282 }
283
284 } // namespace geos.operation.buffer
285 } // namespace geos.operation
286 } // namespace geos
287
288