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