/src/gdal/ogr/ogrcircularstring.cpp
Line | Count | Source (jump to first uncovered line) |
1 | | /****************************************************************************** |
2 | | * |
3 | | * Project: OpenGIS Simple Features Reference Implementation |
4 | | * Purpose: The OGRCircularString geometry class. |
5 | | * Author: Even Rouault, even dot rouault at spatialys dot com |
6 | | * |
7 | | ****************************************************************************** |
8 | | * Copyright (c) 2010, 2014, Even Rouault <even dot rouault at spatialys dot |
9 | | *com> |
10 | | * |
11 | | * SPDX-License-Identifier: MIT |
12 | | ****************************************************************************/ |
13 | | |
14 | | #include "cpl_port.h" |
15 | | #include "ogr_geometry.h" |
16 | | |
17 | | #include <cmath> |
18 | | #include <cstring> |
19 | | |
20 | | #include <algorithm> |
21 | | #include <limits> |
22 | | #include <vector> |
23 | | |
24 | | #include "cpl_error.h" |
25 | | #include "ogr_core.h" |
26 | | #include "ogr_geometry.h" |
27 | | #include "ogr_p.h" |
28 | | |
29 | | static inline double dist(double x0, double y0, double x1, double y1) |
30 | 0 | { |
31 | 0 | return std::sqrt((x1 - x0) * (x1 - x0) + (y1 - y0) * (y1 - y0)); |
32 | 0 | } |
33 | | |
34 | | /************************************************************************/ |
35 | | /* OGRCircularString( const OGRCircularString& ) */ |
36 | | /************************************************************************/ |
37 | | |
38 | | /** |
39 | | * \brief Copy constructor. |
40 | | * |
41 | | * Note: before GDAL 2.1, only the default implementation of the constructor |
42 | | * existed, which could be unsafe to use. |
43 | | * |
44 | | * @since GDAL 2.1 |
45 | | */ |
46 | | |
47 | 0 | OGRCircularString::OGRCircularString(const OGRCircularString &) = default; |
48 | | |
49 | | /************************************************************************/ |
50 | | /* operator=( const OGRCircularString& ) */ |
51 | | /************************************************************************/ |
52 | | |
53 | | /** |
54 | | * \brief Assignment operator. |
55 | | * |
56 | | * Note: before GDAL 2.1, only the default implementation of the operator |
57 | | * existed, which could be unsafe to use. |
58 | | * |
59 | | * @since GDAL 2.1 |
60 | | */ |
61 | | |
62 | | OGRCircularString &OGRCircularString::operator=(const OGRCircularString &other) |
63 | 0 | { |
64 | 0 | if (this != &other) |
65 | 0 | { |
66 | 0 | OGRSimpleCurve::operator=(other); |
67 | 0 | } |
68 | 0 | return *this; |
69 | 0 | } |
70 | | |
71 | | /************************************************************************/ |
72 | | /* getGeometryType() */ |
73 | | /************************************************************************/ |
74 | | |
75 | | OGRwkbGeometryType OGRCircularString::getGeometryType() const |
76 | | |
77 | 12 | { |
78 | 12 | if ((flags & OGR_G_3D) && (flags & OGR_G_MEASURED)) |
79 | 6 | return wkbCircularStringZM; |
80 | 6 | else if (flags & OGR_G_MEASURED) |
81 | 0 | return wkbCircularStringM; |
82 | 6 | else if (flags & OGR_G_3D) |
83 | 5 | return wkbCircularStringZ; |
84 | 1 | else |
85 | 1 | return wkbCircularString; |
86 | 12 | } |
87 | | |
88 | | /************************************************************************/ |
89 | | /* getGeometryName() */ |
90 | | /************************************************************************/ |
91 | | |
92 | | const char *OGRCircularString::getGeometryName() const |
93 | | |
94 | 18 | { |
95 | 18 | return "CIRCULARSTRING"; |
96 | 18 | } |
97 | | |
98 | | /************************************************************************/ |
99 | | /* clone() */ |
100 | | /************************************************************************/ |
101 | | |
102 | | OGRCircularString *OGRCircularString::clone() const |
103 | | |
104 | 0 | { |
105 | 0 | return new (std::nothrow) OGRCircularString(*this); |
106 | 0 | } |
107 | | |
108 | | /************************************************************************/ |
109 | | /* importFromWkb() */ |
110 | | /* */ |
111 | | /* Initialize from serialized stream in well known binary */ |
112 | | /* format. */ |
113 | | /************************************************************************/ |
114 | | |
115 | | OGRErr OGRCircularString::importFromWkb(const unsigned char *pabyData, |
116 | | size_t nSize, OGRwkbVariant eWkbVariant, |
117 | | size_t &nBytesConsumedOut) |
118 | | |
119 | 0 | { |
120 | 0 | OGRErr eErr = OGRSimpleCurve::importFromWkb(pabyData, nSize, eWkbVariant, |
121 | 0 | nBytesConsumedOut); |
122 | 0 | if (eErr == OGRERR_NONE) |
123 | 0 | { |
124 | 0 | if (!IsValidFast()) |
125 | 0 | { |
126 | 0 | empty(); |
127 | 0 | return OGRERR_CORRUPT_DATA; |
128 | 0 | } |
129 | 0 | } |
130 | 0 | return eErr; |
131 | 0 | } |
132 | | |
133 | | /************************************************************************/ |
134 | | /* exportToWkb() */ |
135 | | /* */ |
136 | | /* Build a well known binary representation of this object. */ |
137 | | /************************************************************************/ |
138 | | |
139 | | OGRErr |
140 | | OGRCircularString::exportToWkb(unsigned char *pabyData, |
141 | | const OGRwkbExportOptions *psOptions) const |
142 | | |
143 | 0 | { |
144 | 0 | if (!IsValidFast()) |
145 | 0 | { |
146 | 0 | return OGRERR_FAILURE; |
147 | 0 | } |
148 | | |
149 | 0 | OGRwkbExportOptions sOptions(psOptions ? *psOptions |
150 | 0 | : OGRwkbExportOptions()); |
151 | | |
152 | | // Does not make sense for new geometries, so patch it. |
153 | 0 | if (sOptions.eWkbVariant == wkbVariantOldOgc) |
154 | 0 | sOptions.eWkbVariant = wkbVariantIso; |
155 | 0 | return OGRSimpleCurve::exportToWkb(pabyData, &sOptions); |
156 | 0 | } |
157 | | |
158 | | /************************************************************************/ |
159 | | /* importFromWkt() */ |
160 | | /* */ |
161 | | /* Instantiate from well known text format. Currently this is */ |
162 | | /* `CIRCULARSTRING [Z] ( x y [z], x y [z], ...)', */ |
163 | | /************************************************************************/ |
164 | | |
165 | | OGRErr OGRCircularString::importFromWkt(const char **ppszInput) |
166 | | |
167 | 18 | { |
168 | 18 | const OGRErr eErr = OGRSimpleCurve::importFromWkt(ppszInput); |
169 | 18 | if (eErr == OGRERR_NONE) |
170 | 14 | { |
171 | 14 | if (!IsValidFast()) |
172 | 2 | { |
173 | 2 | empty(); |
174 | 2 | return OGRERR_CORRUPT_DATA; |
175 | 2 | } |
176 | 14 | } |
177 | 16 | return eErr; |
178 | 18 | } |
179 | | |
180 | | /************************************************************************/ |
181 | | /* exportToWkt() */ |
182 | | /************************************************************************/ |
183 | | |
184 | | std::string OGRCircularString::exportToWkt(const OGRWktOptions &opts, |
185 | | OGRErr *err) const |
186 | 0 | { |
187 | 0 | if (!IsValidFast()) |
188 | 0 | { |
189 | 0 | if (err) |
190 | 0 | *err = OGRERR_FAILURE; |
191 | 0 | return std::string(); |
192 | 0 | } |
193 | | |
194 | 0 | OGRWktOptions optsModified(opts); |
195 | 0 | optsModified.variant = wkbVariantIso; |
196 | 0 | return OGRSimpleCurve::exportToWkt(optsModified, err); |
197 | 0 | } |
198 | | |
199 | | /************************************************************************/ |
200 | | /* get_Length() */ |
201 | | /* */ |
202 | | /* For now we return a simple euclidean 2D distance. */ |
203 | | /************************************************************************/ |
204 | | |
205 | | double OGRCircularString::get_Length() const |
206 | | |
207 | 0 | { |
208 | 0 | double dfLength = 0.0; |
209 | 0 | for (int i = 0; i < nPointCount - 2; i += 2) |
210 | 0 | { |
211 | 0 | const double x0 = paoPoints[i].x; |
212 | 0 | const double y0 = paoPoints[i].y; |
213 | 0 | const double x1 = paoPoints[i + 1].x; |
214 | 0 | const double y1 = paoPoints[i + 1].y; |
215 | 0 | const double x2 = paoPoints[i + 2].x; |
216 | 0 | const double y2 = paoPoints[i + 2].y; |
217 | 0 | double R = 0.0; |
218 | 0 | double cx = 0.0; |
219 | 0 | double cy = 0.0; |
220 | 0 | double alpha0 = 0.0; |
221 | 0 | double alpha1 = 0.0; |
222 | 0 | double alpha2 = 0.0; |
223 | 0 | if (OGRGeometryFactory::GetCurveParameters( |
224 | 0 | x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2)) |
225 | 0 | { |
226 | 0 | dfLength += fabs(alpha2 - alpha0) * R; |
227 | 0 | } |
228 | 0 | else |
229 | 0 | { |
230 | 0 | dfLength += dist(x0, y0, x2, y2); |
231 | 0 | } |
232 | 0 | } |
233 | 0 | return dfLength; |
234 | 0 | } |
235 | | |
236 | | /************************************************************************/ |
237 | | /* ExtendEnvelopeWithCircular() */ |
238 | | /************************************************************************/ |
239 | | |
240 | | void OGRCircularString::ExtendEnvelopeWithCircular( |
241 | | OGREnvelope *psEnvelope) const |
242 | 0 | { |
243 | 0 | if (!IsValidFast() || nPointCount == 0) |
244 | 0 | return; |
245 | | |
246 | | // Loop through circular portions and determine if they include some |
247 | | // extremities of the circle. |
248 | 0 | for (int i = 0; i < nPointCount - 2; i += 2) |
249 | 0 | { |
250 | 0 | const double x0 = paoPoints[i].x; |
251 | 0 | const double y0 = paoPoints[i].y; |
252 | 0 | const double x1 = paoPoints[i + 1].x; |
253 | 0 | const double y1 = paoPoints[i + 1].y; |
254 | 0 | const double x2 = paoPoints[i + 2].x; |
255 | 0 | const double y2 = paoPoints[i + 2].y; |
256 | 0 | double R = 0.0; |
257 | 0 | double cx = 0.0; |
258 | 0 | double cy = 0.0; |
259 | 0 | double alpha0 = 0.0; |
260 | 0 | double alpha1 = 0.0; |
261 | 0 | double alpha2 = 0.0; |
262 | 0 | if (OGRGeometryFactory::GetCurveParameters( |
263 | 0 | x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2)) |
264 | 0 | { |
265 | 0 | if (std::isnan(alpha0) || std::isnan(alpha2)) |
266 | 0 | { |
267 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
268 | 0 | "GetCurveParameters returned NaN"); |
269 | 0 | continue; |
270 | 0 | } |
271 | 0 | int quadrantStart = |
272 | 0 | static_cast<int>(std::floor(alpha0 / (M_PI / 2))); |
273 | 0 | int quadrantEnd = static_cast<int>(std::floor(alpha2 / (M_PI / 2))); |
274 | 0 | if (quadrantStart > quadrantEnd) |
275 | 0 | { |
276 | 0 | std::swap(quadrantStart, quadrantEnd); |
277 | 0 | } |
278 | | // Transition trough quadrants in counter-clock wise direction. |
279 | 0 | for (int j = quadrantStart + 1; j <= quadrantEnd; ++j) |
280 | 0 | { |
281 | 0 | switch ((j + 8) % 4) |
282 | 0 | { |
283 | 0 | case 0: |
284 | 0 | psEnvelope->MaxX = std::max(psEnvelope->MaxX, cx + R); |
285 | 0 | break; |
286 | 0 | case 1: |
287 | 0 | psEnvelope->MaxY = std::max(psEnvelope->MaxY, cy + R); |
288 | 0 | break; |
289 | 0 | case 2: |
290 | 0 | psEnvelope->MinX = std::min(psEnvelope->MinX, cx - R); |
291 | 0 | break; |
292 | 0 | case 3: |
293 | 0 | psEnvelope->MinY = std::min(psEnvelope->MaxY, cy - R); |
294 | 0 | break; |
295 | 0 | default: |
296 | 0 | CPLAssert(false); |
297 | 0 | break; |
298 | 0 | } |
299 | 0 | } |
300 | 0 | } |
301 | 0 | } |
302 | 0 | } |
303 | | |
304 | | /************************************************************************/ |
305 | | /* getEnvelope() */ |
306 | | /************************************************************************/ |
307 | | |
308 | | void OGRCircularString::getEnvelope(OGREnvelope *psEnvelope) const |
309 | | |
310 | 0 | { |
311 | 0 | OGRSimpleCurve::getEnvelope(psEnvelope); |
312 | 0 | ExtendEnvelopeWithCircular(psEnvelope); |
313 | 0 | } |
314 | | |
315 | | /************************************************************************/ |
316 | | /* getEnvelope() */ |
317 | | /************************************************************************/ |
318 | | |
319 | | void OGRCircularString::getEnvelope(OGREnvelope3D *psEnvelope) const |
320 | | |
321 | 0 | { |
322 | 0 | OGRSimpleCurve::getEnvelope(psEnvelope); |
323 | 0 | ExtendEnvelopeWithCircular(psEnvelope); |
324 | 0 | } |
325 | | |
326 | | /************************************************************************/ |
327 | | /* OGRCircularString::segmentize() */ |
328 | | /************************************************************************/ |
329 | | |
330 | | bool OGRCircularString::segmentize(double dfMaxLength) |
331 | 0 | { |
332 | 0 | if (!IsValidFast()) |
333 | 0 | return false; |
334 | 0 | if (nPointCount == 0) |
335 | 0 | return true; |
336 | | |
337 | | // So as to make sure that the same line followed in both directions |
338 | | // result in the same segmentized line. |
339 | 0 | if (paoPoints[0].x < paoPoints[nPointCount - 1].x || |
340 | 0 | (paoPoints[0].x == paoPoints[nPointCount - 1].x && |
341 | 0 | paoPoints[0].y < paoPoints[nPointCount - 1].y)) |
342 | 0 | { |
343 | 0 | reversePoints(); |
344 | 0 | const bool bRet = segmentize(dfMaxLength); |
345 | 0 | reversePoints(); |
346 | 0 | return bRet; |
347 | 0 | } |
348 | | |
349 | 0 | std::vector<OGRRawPoint> aoRawPoint; |
350 | 0 | std::vector<double> adfZ; |
351 | 0 | bool bRet = true; |
352 | 0 | for (int i = 0; i < nPointCount - 2; i += 2) |
353 | 0 | { |
354 | 0 | const double x0 = paoPoints[i].x; |
355 | 0 | const double y0 = paoPoints[i].y; |
356 | 0 | const double x1 = paoPoints[i + 1].x; |
357 | 0 | const double y1 = paoPoints[i + 1].y; |
358 | 0 | const double x2 = paoPoints[i + 2].x; |
359 | 0 | const double y2 = paoPoints[i + 2].y; |
360 | 0 | double R = 0.0; |
361 | 0 | double cx = 0.0; |
362 | 0 | double cy = 0.0; |
363 | 0 | double alpha0 = 0.0; |
364 | 0 | double alpha1 = 0.0; |
365 | 0 | double alpha2 = 0.0; |
366 | |
|
367 | 0 | aoRawPoint.emplace_back(x0, y0); |
368 | 0 | if (padfZ) |
369 | 0 | adfZ.emplace_back(padfZ[i]); |
370 | | |
371 | | // We have strong constraints on the number of intermediate points |
372 | | // we can add. |
373 | |
|
374 | 0 | if (OGRGeometryFactory::GetCurveParameters( |
375 | 0 | x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2)) |
376 | 0 | { |
377 | | // It is an arc circle. |
378 | 0 | const double dfSegmentLength1 = fabs(alpha1 - alpha0) * R; |
379 | 0 | const double dfSegmentLength2 = fabs(alpha2 - alpha1) * R; |
380 | 0 | if (dfSegmentLength1 > dfMaxLength || |
381 | 0 | dfSegmentLength2 > dfMaxLength) |
382 | 0 | { |
383 | 0 | const double dfVal = |
384 | 0 | 1 + 2 * std::floor(dfSegmentLength1 / dfMaxLength / 2.0); |
385 | 0 | if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 || |
386 | 0 | std::isnan(dfVal)) |
387 | 0 | { |
388 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
389 | 0 | "segmentize nIntermediatePoints invalid: %lf", |
390 | 0 | dfVal); |
391 | 0 | bRet = false; |
392 | 0 | break; |
393 | 0 | } |
394 | 0 | const int nIntermediatePoints = static_cast<int>(dfVal); |
395 | 0 | const double dfStep = |
396 | 0 | (alpha1 - alpha0) / (nIntermediatePoints + 1); |
397 | 0 | for (int j = 1; j <= nIntermediatePoints; ++j) |
398 | 0 | { |
399 | 0 | double alpha = alpha0 + dfStep * j; |
400 | 0 | const double x = cx + R * cos(alpha); |
401 | 0 | const double y = cy + R * sin(alpha); |
402 | 0 | aoRawPoint.emplace_back(x, y); |
403 | 0 | if (padfZ) |
404 | 0 | { |
405 | 0 | const double z = padfZ[i] + (padfZ[i + 1] - padfZ[i]) * |
406 | 0 | (alpha - alpha0) / |
407 | 0 | (alpha1 - alpha0); |
408 | 0 | adfZ.emplace_back(z); |
409 | 0 | } |
410 | 0 | } |
411 | 0 | } |
412 | 0 | aoRawPoint.emplace_back(x1, y1); |
413 | 0 | if (padfZ) |
414 | 0 | adfZ.emplace_back(padfZ[i + 1]); |
415 | |
|
416 | 0 | if (dfSegmentLength1 > dfMaxLength || |
417 | 0 | dfSegmentLength2 > dfMaxLength) |
418 | 0 | { |
419 | 0 | const double dfVal = |
420 | 0 | 1 + 2 * std::floor(dfSegmentLength2 / dfMaxLength / 2.0); |
421 | 0 | if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 || |
422 | 0 | std::isnan(dfVal)) |
423 | 0 | { |
424 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
425 | 0 | "segmentize nIntermediatePoints invalid 2: %lf", |
426 | 0 | dfVal); |
427 | 0 | bRet = false; |
428 | 0 | break; |
429 | 0 | } |
430 | 0 | int nIntermediatePoints = static_cast<int>(dfVal); |
431 | 0 | const double dfStep = |
432 | 0 | (alpha2 - alpha1) / (nIntermediatePoints + 1); |
433 | 0 | for (int j = 1; j <= nIntermediatePoints; ++j) |
434 | 0 | { |
435 | 0 | const double alpha = alpha1 + dfStep * j; |
436 | 0 | const double x = cx + R * cos(alpha); |
437 | 0 | const double y = cy + R * sin(alpha); |
438 | 0 | aoRawPoint.emplace_back(x, y); |
439 | 0 | if (padfZ) |
440 | 0 | { |
441 | 0 | const double z = |
442 | 0 | padfZ[i + 1] + (padfZ[i + 2] - padfZ[i + 1]) * |
443 | 0 | (alpha - alpha1) / |
444 | 0 | (alpha2 - alpha1); |
445 | 0 | adfZ.emplace_back(z); |
446 | 0 | } |
447 | 0 | } |
448 | 0 | } |
449 | 0 | } |
450 | 0 | else |
451 | 0 | { |
452 | | // It is a straight line. |
453 | 0 | const double dfSegmentLength1 = dist(x0, y0, x1, y1); |
454 | 0 | const double dfSegmentLength2 = dist(x1, y1, x2, y2); |
455 | 0 | if (dfSegmentLength1 > dfMaxLength || |
456 | 0 | dfSegmentLength2 > dfMaxLength) |
457 | 0 | { |
458 | 0 | const double dfVal = |
459 | 0 | 1 + 2 * std::ceil(dfSegmentLength1 / dfMaxLength / 2.0); |
460 | 0 | if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 || |
461 | 0 | std::isnan(dfVal)) |
462 | 0 | { |
463 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
464 | 0 | "segmentize nIntermediatePoints invalid 2: %lf", |
465 | 0 | dfVal); |
466 | 0 | bRet = false; |
467 | 0 | break; |
468 | 0 | } |
469 | 0 | int nIntermediatePoints = static_cast<int>(dfVal); |
470 | 0 | for (int j = 1; j <= nIntermediatePoints; ++j) |
471 | 0 | { |
472 | 0 | aoRawPoint.emplace_back( |
473 | 0 | x0 + j * (x1 - x0) / (nIntermediatePoints + 1), |
474 | 0 | y0 + j * (y1 - y0) / (nIntermediatePoints + 1)); |
475 | 0 | if (padfZ) |
476 | 0 | adfZ.emplace_back(padfZ[i] + |
477 | 0 | j * (padfZ[i + 1] - padfZ[i]) / |
478 | 0 | (nIntermediatePoints + 1)); |
479 | 0 | } |
480 | 0 | } |
481 | | |
482 | 0 | aoRawPoint.emplace_back(x1, y1); |
483 | 0 | if (padfZ) |
484 | 0 | adfZ.emplace_back(padfZ[i + 1]); |
485 | |
|
486 | 0 | if (dfSegmentLength1 > dfMaxLength || |
487 | 0 | dfSegmentLength2 > dfMaxLength) |
488 | 0 | { |
489 | 0 | const double dfVal = |
490 | 0 | 1 + 2 * std::ceil(dfSegmentLength2 / dfMaxLength / 2.0); |
491 | 0 | if (dfVal >= std::numeric_limits<int>::max() || dfVal < 0.0 || |
492 | 0 | std::isnan(dfVal)) |
493 | 0 | { |
494 | 0 | CPLError(CE_Failure, CPLE_AppDefined, |
495 | 0 | "segmentize nIntermediatePoints invalid 3: %lf", |
496 | 0 | dfVal); |
497 | 0 | bRet = false; |
498 | 0 | break; |
499 | 0 | } |
500 | 0 | const int nIntermediatePoints = static_cast<int>(dfVal); |
501 | |
|
502 | 0 | for (int j = 1; j <= nIntermediatePoints; ++j) |
503 | 0 | { |
504 | 0 | aoRawPoint.emplace_back( |
505 | 0 | x1 + j * (x2 - x1) / (nIntermediatePoints + 1), |
506 | 0 | y1 + j * (y2 - y1) / (nIntermediatePoints + 1)); |
507 | 0 | if (padfZ) |
508 | 0 | adfZ.emplace_back(padfZ[i + 1] + |
509 | 0 | j * (padfZ[i + 2] - padfZ[i + 1]) / |
510 | 0 | (nIntermediatePoints + 1)); |
511 | 0 | } |
512 | 0 | } |
513 | 0 | } |
514 | 0 | } |
515 | 0 | aoRawPoint.push_back(paoPoints[nPointCount - 1]); |
516 | 0 | if (padfZ) |
517 | 0 | adfZ.push_back(padfZ[nPointCount - 1]); |
518 | |
|
519 | 0 | CPLAssert(aoRawPoint.empty() || |
520 | 0 | (aoRawPoint.size() >= 3 && (aoRawPoint.size() % 2) == 1)); |
521 | 0 | if (padfZ) |
522 | 0 | { |
523 | 0 | CPLAssert(adfZ.size() == aoRawPoint.size()); |
524 | 0 | } |
525 | | |
526 | | // Is there actually something to modify? |
527 | 0 | if (nPointCount < static_cast<int>(aoRawPoint.size())) |
528 | 0 | { |
529 | 0 | nPointCount = static_cast<int>(aoRawPoint.size()); |
530 | 0 | paoPoints = static_cast<OGRRawPoint *>( |
531 | 0 | CPLRealloc(paoPoints, sizeof(OGRRawPoint) * nPointCount)); |
532 | 0 | memcpy(paoPoints, &aoRawPoint[0], sizeof(OGRRawPoint) * nPointCount); |
533 | 0 | if (padfZ) |
534 | 0 | { |
535 | 0 | padfZ = static_cast<double *>( |
536 | 0 | CPLRealloc(padfZ, sizeof(double) * aoRawPoint.size())); |
537 | 0 | memcpy(padfZ, &adfZ[0], sizeof(double) * nPointCount); |
538 | 0 | } |
539 | 0 | } |
540 | 0 | return bRet; |
541 | 0 | } |
542 | | |
543 | | /************************************************************************/ |
544 | | /* Value() */ |
545 | | /* */ |
546 | | /* Get an interpolated point at some distance along the curve. */ |
547 | | /************************************************************************/ |
548 | | |
549 | | void OGRCircularString::Value(double dfDistance, OGRPoint *poPoint) const |
550 | | |
551 | 0 | { |
552 | 0 | if (dfDistance < 0) |
553 | 0 | { |
554 | 0 | StartPoint(poPoint); |
555 | 0 | return; |
556 | 0 | } |
557 | | |
558 | 0 | double dfLength = 0; |
559 | |
|
560 | 0 | for (int i = 0; i < nPointCount - 2; i += 2) |
561 | 0 | { |
562 | 0 | const double x0 = paoPoints[i].x; |
563 | 0 | const double y0 = paoPoints[i].y; |
564 | 0 | const double x1 = paoPoints[i + 1].x; |
565 | 0 | const double y1 = paoPoints[i + 1].y; |
566 | 0 | const double x2 = paoPoints[i + 2].x; |
567 | 0 | const double y2 = paoPoints[i + 2].y; |
568 | 0 | double R = 0.0; |
569 | 0 | double cx = 0.0; |
570 | 0 | double cy = 0.0; |
571 | 0 | double alpha0 = 0.0; |
572 | 0 | double alpha1 = 0.0; |
573 | 0 | double alpha2 = 0.0; |
574 | | |
575 | | // We have strong constraints on the number of intermediate points |
576 | | // we can add. |
577 | |
|
578 | 0 | if (OGRGeometryFactory::GetCurveParameters( |
579 | 0 | x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2)) |
580 | 0 | { |
581 | | // It is an arc circle. |
582 | 0 | const double dfSegLength = fabs(alpha2 - alpha0) * R; |
583 | 0 | if (dfSegLength > 0) |
584 | 0 | { |
585 | 0 | if ((dfLength <= dfDistance) && |
586 | 0 | ((dfLength + dfSegLength) >= dfDistance)) |
587 | 0 | { |
588 | 0 | const double dfRatio = |
589 | 0 | (dfDistance - dfLength) / dfSegLength; |
590 | |
|
591 | 0 | const double alpha = |
592 | 0 | alpha0 * (1 - dfRatio) + alpha2 * dfRatio; |
593 | 0 | const double x = cx + R * cos(alpha); |
594 | 0 | const double y = cy + R * sin(alpha); |
595 | |
|
596 | 0 | poPoint->setX(x); |
597 | 0 | poPoint->setY(y); |
598 | |
|
599 | 0 | if (getCoordinateDimension() == 3) |
600 | 0 | poPoint->setZ(padfZ[i] * (1 - dfRatio) + |
601 | 0 | padfZ[i + 2] * dfRatio); |
602 | |
|
603 | 0 | return; |
604 | 0 | } |
605 | | |
606 | 0 | dfLength += dfSegLength; |
607 | 0 | } |
608 | 0 | } |
609 | 0 | else |
610 | 0 | { |
611 | | // It is a straight line. |
612 | 0 | const double dfSegLength = dist(x0, y0, x2, y2); |
613 | 0 | if (dfSegLength > 0) |
614 | 0 | { |
615 | 0 | if ((dfLength <= dfDistance) && |
616 | 0 | ((dfLength + dfSegLength) >= dfDistance)) |
617 | 0 | { |
618 | 0 | const double dfRatio = |
619 | 0 | (dfDistance - dfLength) / dfSegLength; |
620 | |
|
621 | 0 | poPoint->setX(paoPoints[i].x * (1 - dfRatio) + |
622 | 0 | paoPoints[i + 2].x * dfRatio); |
623 | 0 | poPoint->setY(paoPoints[i].y * (1 - dfRatio) + |
624 | 0 | paoPoints[i + 2].y * dfRatio); |
625 | |
|
626 | 0 | if (getCoordinateDimension() == 3) |
627 | 0 | poPoint->setZ(padfZ[i] * (1 - dfRatio) + |
628 | 0 | padfZ[i + 2] * dfRatio); |
629 | |
|
630 | 0 | return; |
631 | 0 | } |
632 | | |
633 | 0 | dfLength += dfSegLength; |
634 | 0 | } |
635 | 0 | } |
636 | 0 | } |
637 | | |
638 | 0 | EndPoint(poPoint); |
639 | 0 | } |
640 | | |
641 | | /************************************************************************/ |
642 | | /* CurveToLine() */ |
643 | | /************************************************************************/ |
644 | | |
645 | | OGRLineString * |
646 | | OGRCircularString::CurveToLine(double dfMaxAngleStepSizeDegrees, |
647 | | const char *const *papszOptions) const |
648 | 0 | { |
649 | 0 | OGRLineString *poLine = new OGRLineString(); |
650 | 0 | poLine->assignSpatialReference(getSpatialReference()); |
651 | |
|
652 | 0 | const bool bHasZ = getCoordinateDimension() == 3; |
653 | 0 | for (int i = 0; i < nPointCount - 2; i += 2) |
654 | 0 | { |
655 | 0 | OGRLineString *poArc = OGRGeometryFactory::curveToLineString( |
656 | 0 | paoPoints[i].x, paoPoints[i].y, padfZ ? padfZ[i] : 0.0, |
657 | 0 | paoPoints[i + 1].x, paoPoints[i + 1].y, padfZ ? padfZ[i + 1] : 0.0, |
658 | 0 | paoPoints[i + 2].x, paoPoints[i + 2].y, padfZ ? padfZ[i + 2] : 0.0, |
659 | 0 | bHasZ, dfMaxAngleStepSizeDegrees, papszOptions); |
660 | 0 | poLine->addSubLineString(poArc, (i == 0) ? 0 : 1); |
661 | 0 | delete poArc; |
662 | 0 | } |
663 | 0 | return poLine; |
664 | 0 | } |
665 | | |
666 | | /************************************************************************/ |
667 | | /* IsValidFast() */ |
668 | | /************************************************************************/ |
669 | | |
670 | | OGRBoolean OGRCircularString::IsValidFast() const |
671 | | |
672 | 14 | { |
673 | 14 | if (nPointCount == 1 || nPointCount == 2 || |
674 | 14 | (nPointCount >= 3 && (nPointCount % 2) == 0)) |
675 | 2 | { |
676 | 2 | CPLError(CE_Failure, CPLE_NotSupported, |
677 | 2 | "Bad number of points in circular string : %d", nPointCount); |
678 | 2 | return FALSE; |
679 | 2 | } |
680 | 12 | return TRUE; |
681 | 14 | } |
682 | | |
683 | | /************************************************************************/ |
684 | | /* IsValid() */ |
685 | | /************************************************************************/ |
686 | | |
687 | | OGRBoolean OGRCircularString::IsValid() const |
688 | | |
689 | 0 | { |
690 | 0 | return IsValidFast() && OGRGeometry::IsValid(); |
691 | 0 | } |
692 | | |
693 | | /************************************************************************/ |
694 | | /* hasCurveGeometry() */ |
695 | | /************************************************************************/ |
696 | | |
697 | | OGRBoolean |
698 | | OGRCircularString::hasCurveGeometry(int /* bLookForNonLinear */) const |
699 | 12 | { |
700 | 12 | return TRUE; |
701 | 12 | } |
702 | | |
703 | | /************************************************************************/ |
704 | | /* getLinearGeometry() */ |
705 | | /************************************************************************/ |
706 | | |
707 | | OGRGeometry * |
708 | | OGRCircularString::getLinearGeometry(double dfMaxAngleStepSizeDegrees, |
709 | | const char *const *papszOptions) const |
710 | 0 | { |
711 | 0 | return CurveToLine(dfMaxAngleStepSizeDegrees, papszOptions); |
712 | 0 | } |
713 | | |
714 | | //! @cond Doxygen_Suppress |
715 | | /************************************************************************/ |
716 | | /* GetCasterToLineString() */ |
717 | | /************************************************************************/ |
718 | | |
719 | | static OGRLineString *CasterToLineString(OGRCurve *poGeom) |
720 | 0 | { |
721 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible", |
722 | 0 | poGeom->getGeometryName()); |
723 | 0 | delete poGeom; |
724 | 0 | return nullptr; |
725 | 0 | } |
726 | | |
727 | | OGRCurveCasterToLineString OGRCircularString::GetCasterToLineString() const |
728 | 0 | { |
729 | 0 | return ::CasterToLineString; |
730 | 0 | } |
731 | | |
732 | | /************************************************************************/ |
733 | | /* GetCasterToLinearRing() */ |
734 | | /************************************************************************/ |
735 | | |
736 | | static OGRLinearRing *CasterToLinearRing(OGRCurve *poGeom) |
737 | 0 | { |
738 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible", |
739 | 0 | poGeom->getGeometryName()); |
740 | 0 | delete poGeom; |
741 | 0 | return nullptr; |
742 | 0 | } |
743 | | |
744 | | OGRCurveCasterToLinearRing OGRCircularString::GetCasterToLinearRing() const |
745 | 0 | { |
746 | 0 | return ::CasterToLinearRing; |
747 | 0 | } |
748 | | |
749 | | //! @endcond |
750 | | |
751 | | /************************************************************************/ |
752 | | /* IsFullCircle() */ |
753 | | /************************************************************************/ |
754 | | |
755 | | int OGRCircularString::IsFullCircle(double &cx, double &cy, |
756 | | double &square_R) const |
757 | 0 | { |
758 | 0 | if (getNumPoints() == 3 && get_IsClosed()) |
759 | 0 | { |
760 | 0 | const double x0 = getX(0); |
761 | 0 | const double y0 = getY(0); |
762 | 0 | const double x1 = getX(1); |
763 | 0 | const double y1 = getY(1); |
764 | 0 | cx = (x0 + x1) / 2; |
765 | 0 | cy = (y0 + y1) / 2; |
766 | 0 | square_R = (x1 - cx) * (x1 - cx) + (y1 - cy) * (y1 - cy); |
767 | 0 | return TRUE; |
768 | 0 | } |
769 | | // Full circle defined by 2 arcs? |
770 | 0 | else if (getNumPoints() == 5 && get_IsClosed()) |
771 | 0 | { |
772 | 0 | double R_1 = 0.0; |
773 | 0 | double cx_1 = 0.0; |
774 | 0 | double cy_1 = 0.0; |
775 | 0 | double alpha0_1 = 0.0; |
776 | 0 | double alpha1_1 = 0.0; |
777 | 0 | double alpha2_1 = 0.0; |
778 | 0 | double R_2 = 0.0; |
779 | 0 | double cx_2 = 0.0; |
780 | 0 | double cy_2 = 0.0; |
781 | 0 | double alpha0_2 = 0.0; |
782 | 0 | double alpha1_2 = 0.0; |
783 | 0 | double alpha2_2 = 0.0; |
784 | 0 | if (OGRGeometryFactory::GetCurveParameters( |
785 | 0 | getX(0), getY(0), getX(1), getY(1), getX(2), getY(2), R_1, cx_1, |
786 | 0 | cy_1, alpha0_1, alpha1_1, alpha2_1) && |
787 | 0 | OGRGeometryFactory::GetCurveParameters( |
788 | 0 | getX(2), getY(2), getX(3), getY(3), getX(4), getY(4), R_2, cx_2, |
789 | 0 | cy_2, alpha0_2, alpha1_2, alpha2_2) && |
790 | 0 | fabs(R_1 - R_2) < 1e-10 && fabs(cx_1 - cx_2) < 1e-10 && |
791 | 0 | fabs(cy_1 - cy_2) < 1e-10 && |
792 | 0 | (alpha2_1 - alpha0_1) * (alpha2_2 - alpha0_2) > 0) |
793 | 0 | { |
794 | 0 | cx = cx_1; |
795 | 0 | cy = cy_1; |
796 | 0 | square_R = R_1 * R_1; |
797 | 0 | return TRUE; |
798 | 0 | } |
799 | 0 | } |
800 | 0 | return FALSE; |
801 | 0 | } |
802 | | |
803 | | /************************************************************************/ |
804 | | /* get_AreaOfCurveSegments() */ |
805 | | /************************************************************************/ |
806 | | |
807 | | //! @cond Doxygen_Suppress |
808 | | double OGRCircularString::get_AreaOfCurveSegments() const |
809 | 0 | { |
810 | 0 | double dfArea = 0.0; |
811 | 0 | for (int i = 0; i < getNumPoints() - 2; i += 2) |
812 | 0 | { |
813 | 0 | const double x0 = getX(i); |
814 | 0 | const double y0 = getY(i); |
815 | 0 | const double x1 = getX(i + 1); |
816 | 0 | const double y1 = getY(i + 1); |
817 | 0 | const double x2 = getX(i + 2); |
818 | 0 | const double y2 = getY(i + 2); |
819 | 0 | double R = 0.0; |
820 | 0 | double cx = 0.0; |
821 | 0 | double cy = 0.0; |
822 | 0 | double alpha0 = 0.0; |
823 | 0 | double alpha1 = 0.0; |
824 | 0 | double alpha2 = 0.0; |
825 | 0 | if (OGRGeometryFactory::GetCurveParameters( |
826 | 0 | x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2)) |
827 | 0 | { |
828 | | // Should be <= PI in absolute value. |
829 | 0 | const double delta_alpha01 = alpha1 - alpha0; |
830 | 0 | const double delta_alpha12 = alpha2 - alpha1; // Same. |
831 | | // http://en.wikipedia.org/wiki/Circular_segment |
832 | 0 | dfArea += 0.5 * R * R * |
833 | 0 | fabs(delta_alpha01 - sin(delta_alpha01) + delta_alpha12 - |
834 | 0 | sin(delta_alpha12)); |
835 | 0 | } |
836 | 0 | } |
837 | 0 | return dfArea; |
838 | 0 | } |
839 | | |
840 | | //! @endcond |
841 | | |
842 | | /************************************************************************/ |
843 | | /* get_Area() */ |
844 | | /************************************************************************/ |
845 | | |
846 | | double OGRCircularString::get_Area() const |
847 | 0 | { |
848 | 0 | if (IsEmpty() || !get_IsClosed()) |
849 | 0 | return 0; |
850 | | |
851 | 0 | double cx = 0.0; |
852 | 0 | double cy = 0.0; |
853 | 0 | double square_R = 0.0; |
854 | |
|
855 | 0 | if (IsFullCircle(cx, cy, square_R)) |
856 | 0 | { |
857 | 0 | return M_PI * square_R; |
858 | 0 | } |
859 | | |
860 | | // Optimization for convex rings. |
861 | 0 | if (IsConvex()) |
862 | 0 | { |
863 | | // Compute area of shape without the circular segments. |
864 | 0 | double dfArea = get_LinearArea(); |
865 | | |
866 | | // Add the area of the spherical segments. |
867 | 0 | dfArea += get_AreaOfCurveSegments(); |
868 | |
|
869 | 0 | return dfArea; |
870 | 0 | } |
871 | | |
872 | 0 | OGRLineString *poLS = CurveToLine(); |
873 | 0 | const double dfArea = poLS->get_Area(); |
874 | 0 | delete poLS; |
875 | |
|
876 | 0 | return dfArea; |
877 | 0 | } |
878 | | |
879 | | /************************************************************************/ |
880 | | /* get_GeodesicArea() */ |
881 | | /************************************************************************/ |
882 | | |
883 | | double OGRCircularString::get_GeodesicArea( |
884 | | const OGRSpatialReference *poSRSOverride) const |
885 | 0 | { |
886 | 0 | if (IsEmpty()) |
887 | 0 | return 0; |
888 | | |
889 | 0 | if (!get_IsClosed()) |
890 | 0 | { |
891 | 0 | CPLError(CE_Failure, CPLE_AppDefined, "Non-closed geometry"); |
892 | 0 | return -1; |
893 | 0 | } |
894 | | |
895 | 0 | if (!poSRSOverride) |
896 | 0 | poSRSOverride = getSpatialReference(); |
897 | |
|
898 | 0 | auto poLS = std::unique_ptr<OGRLineString>(CurveToLine()); |
899 | 0 | return poLS->get_GeodesicArea(poSRSOverride); |
900 | 0 | } |
901 | | |
902 | | /************************************************************************/ |
903 | | /* get_GeodesicLength() */ |
904 | | /************************************************************************/ |
905 | | |
906 | | double OGRCircularString::get_GeodesicLength( |
907 | | const OGRSpatialReference *poSRSOverride) const |
908 | 0 | { |
909 | 0 | if (IsEmpty()) |
910 | 0 | return 0; |
911 | | |
912 | 0 | if (!poSRSOverride) |
913 | 0 | poSRSOverride = getSpatialReference(); |
914 | |
|
915 | 0 | auto poLS = std::unique_ptr<OGRLineString>(CurveToLine()); |
916 | 0 | return poLS->get_GeodesicLength(poSRSOverride); |
917 | 0 | } |
918 | | |
919 | | //! @cond Doxygen_Suppress |
920 | | |
921 | | /************************************************************************/ |
922 | | /* ContainsPoint() */ |
923 | | /************************************************************************/ |
924 | | |
925 | | int OGRCircularString::ContainsPoint(const OGRPoint *p) const |
926 | 0 | { |
927 | 0 | double cx = 0.0; |
928 | 0 | double cy = 0.0; |
929 | 0 | double square_R = 0.0; |
930 | 0 | if (IsFullCircle(cx, cy, square_R)) |
931 | 0 | { |
932 | 0 | const double square_dist = (p->getX() - cx) * (p->getX() - cx) + |
933 | 0 | (p->getY() - cy) * (p->getY() - cy); |
934 | 0 | return square_dist < square_R; |
935 | 0 | } |
936 | 0 | return -1; |
937 | 0 | } |
938 | | |
939 | | /************************************************************************/ |
940 | | /* IntersectsPoint() */ |
941 | | /************************************************************************/ |
942 | | |
943 | | int OGRCircularString::IntersectsPoint(const OGRPoint *p) const |
944 | 0 | { |
945 | 0 | double cx = 0.0; |
946 | 0 | double cy = 0.0; |
947 | 0 | double square_R = 0.0; |
948 | 0 | if (IsFullCircle(cx, cy, square_R)) |
949 | 0 | { |
950 | 0 | const double square_dist = (p->getX() - cx) * (p->getX() - cx) + |
951 | 0 | (p->getY() - cy) * (p->getY() - cy); |
952 | 0 | return square_dist <= square_R; |
953 | 0 | } |
954 | 0 | return -1; |
955 | 0 | } |
956 | | |
957 | | //! @endcond |