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