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 10586 : OGRwkbGeometryType OGRCircularString::getGeometryType() const
66 :
67 : {
68 10586 : if ((flags & OGR_G_3D) && (flags & OGR_G_MEASURED))
69 683 : return wkbCircularStringZM;
70 9903 : else if (flags & OGR_G_MEASURED)
71 67 : return wkbCircularStringM;
72 9836 : else if (flags & OGR_G_3D)
73 539 : return wkbCircularStringZ;
74 : else
75 9297 : return wkbCircularString;
76 : }
77 :
78 : /************************************************************************/
79 : /* getGeometryName() */
80 : /************************************************************************/
81 :
82 4853 : const char *OGRCircularString::getGeometryName() const
83 :
84 : {
85 4853 : return "CIRCULARSTRING";
86 : }
87 :
88 : /************************************************************************/
89 : /* clone() */
90 : /************************************************************************/
91 :
92 515 : OGRCircularString *OGRCircularString::clone() const
93 :
94 : {
95 515 : 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 3556 : OGRErr OGRCircularString::importFromWkt(const char **ppszInput)
156 :
157 : {
158 3556 : const OGRErr eErr = OGRSimpleCurve::importFromWkt(ppszInput);
159 3556 : if (eErr == OGRERR_NONE)
160 : {
161 3556 : if (!IsValidFast())
162 : {
163 1 : empty();
164 1 : return OGRERR_CORRUPT_DATA;
165 : }
166 : }
167 3555 : 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 257 : void OGRCircularString::ExtendEnvelopeWithCircular(
231 : OGREnvelope *psEnvelope) const
232 : {
233 257 : 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 570 : for (int i = 0; i < nPointCount - 2; i += 2)
239 : {
240 313 : const double x0 = paoPoints[i].x;
241 313 : const double y0 = paoPoints[i].y;
242 313 : const double x1 = paoPoints[i + 1].x;
243 313 : const double y1 = paoPoints[i + 1].y;
244 313 : const double x2 = paoPoints[i + 2].x;
245 313 : const double y2 = paoPoints[i + 2].y;
246 313 : double R = 0.0;
247 313 : double cx = 0.0;
248 313 : double cy = 0.0;
249 313 : double alpha0 = 0.0;
250 313 : double alpha1 = 0.0;
251 313 : double alpha2 = 0.0;
252 313 : if (OGRGeometryFactory::GetCurveParameters(
253 313 : x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
254 : {
255 313 : if (std::isnan(alpha0) || std::isnan(alpha2))
256 : {
257 0 : CPLError(CE_Failure, CPLE_AppDefined,
258 : "GetCurveParameters returned NaN");
259 0 : continue;
260 : }
261 313 : int quadrantStart =
262 313 : static_cast<int>(std::floor(alpha0 / (M_PI / 2)));
263 313 : int quadrantEnd = static_cast<int>(std::floor(alpha2 / (M_PI / 2)));
264 313 : if (quadrantStart > quadrantEnd)
265 : {
266 200 : std::swap(quadrantStart, quadrantEnd);
267 : }
268 : // Transition through quadrants in counter-clock wise direction.
269 1010 : for (int j = quadrantStart + 1; j <= quadrantEnd; ++j)
270 : {
271 697 : switch ((j + 8) % 4)
272 : {
273 153 : case 0:
274 153 : psEnvelope->MaxX = std::max(psEnvelope->MaxX, cx + R);
275 153 : break;
276 221 : case 1:
277 221 : psEnvelope->MaxY = std::max(psEnvelope->MaxY, cy + R);
278 221 : break;
279 203 : case 2:
280 203 : psEnvelope->MinX = std::min(psEnvelope->MinX, cx - R);
281 203 : break;
282 120 : case 3:
283 120 : psEnvelope->MinY = std::min(psEnvelope->MaxY, cy - R);
284 120 : break;
285 0 : default:
286 0 : CPLAssert(false);
287 : break;
288 : }
289 : }
290 : }
291 : }
292 : }
293 :
294 : /************************************************************************/
295 : /* getEnvelope() */
296 : /************************************************************************/
297 :
298 141 : void OGRCircularString::getEnvelope(OGREnvelope *psEnvelope) const
299 :
300 : {
301 141 : OGRSimpleCurve::getEnvelope(psEnvelope);
302 141 : ExtendEnvelopeWithCircular(psEnvelope);
303 141 : }
304 :
305 : /************************************************************************/
306 : /* getEnvelope() */
307 : /************************************************************************/
308 :
309 116 : void OGRCircularString::getEnvelope(OGREnvelope3D *psEnvelope) const
310 :
311 : {
312 116 : OGRSimpleCurve::getEnvelope(psEnvelope);
313 116 : ExtendEnvelopeWithCircular(psEnvelope);
314 116 : }
315 :
316 : /************************************************************************/
317 : /* OGRCircularString::segmentize() */
318 : /************************************************************************/
319 :
320 17 : bool OGRCircularString::segmentize(double dfMaxLength)
321 : {
322 17 : if (!IsValidFast())
323 0 : return false;
324 17 : 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 17 : if (paoPoints[0].x < paoPoints[nPointCount - 1].x ||
330 10 : (paoPoints[0].x == paoPoints[nPointCount - 1].x &&
331 2 : paoPoints[0].y < paoPoints[nPointCount - 1].y))
332 : {
333 7 : reversePoints();
334 7 : const bool bRet = segmentize(dfMaxLength);
335 7 : reversePoints();
336 7 : return bRet;
337 : }
338 :
339 20 : std::vector<OGRRawPoint> aoRawPoint;
340 20 : std::vector<double> adfZ;
341 20 : std::vector<double> adfM;
342 10 : bool bRet = true;
343 21 : for (int i = 0; i < nPointCount - 2; i += 2)
344 : {
345 12 : const double x0 = paoPoints[i].x;
346 12 : const double y0 = paoPoints[i].y;
347 12 : const double x1 = paoPoints[i + 1].x;
348 12 : const double y1 = paoPoints[i + 1].y;
349 12 : const double x2 = paoPoints[i + 2].x;
350 12 : const double y2 = paoPoints[i + 2].y;
351 12 : double R = 0.0;
352 12 : double cx = 0.0;
353 12 : double cy = 0.0;
354 12 : double alpha0 = 0.0;
355 12 : double alpha1 = 0.0;
356 12 : double alpha2 = 0.0;
357 :
358 12 : aoRawPoint.emplace_back(x0, y0);
359 12 : if (padfZ)
360 6 : adfZ.emplace_back(padfZ[i]);
361 12 : if (padfM)
362 4 : adfM.emplace_back(padfM[i]);
363 :
364 12 : constexpr int kMax = 2 << 26;
365 :
366 : // We have strong constraints on the number of intermediate points
367 : // we can add.
368 :
369 12 : if (OGRGeometryFactory::GetCurveParameters(
370 12 : x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
371 : {
372 : // It is an arc circle.
373 10 : const double dfSegmentLength1 = fabs(alpha1 - alpha0) * R;
374 10 : const double dfSegmentLength2 = fabs(alpha2 - alpha1) * R;
375 10 : if (dfSegmentLength1 > dfMaxLength ||
376 : dfSegmentLength2 > dfMaxLength)
377 : {
378 10 : const double dfVal =
379 10 : 1 + 2 * std::floor(dfSegmentLength1 / dfMaxLength / 2.0);
380 20 : if (dfVal < 0.0 || std::isnan(dfVal) ||
381 10 : dfVal >= kMax - static_cast<int>(aoRawPoint.size()))
382 : {
383 1 : CPLError(CE_Failure, CPLE_AppDefined,
384 : "segmentize nIntermediatePoints invalid: %lf",
385 : dfVal);
386 1 : return false;
387 : }
388 9 : const int nIntermediatePoints = static_cast<int>(dfVal);
389 9 : const double dfStep =
390 9 : (alpha1 - alpha0) / (nIntermediatePoints + 1);
391 46 : for (int j = 1; j <= nIntermediatePoints; ++j)
392 : {
393 37 : double alpha = alpha0 + dfStep * j;
394 37 : const double x = cx + R * cos(alpha);
395 37 : const double y = cy + R * sin(alpha);
396 37 : aoRawPoint.emplace_back(x, y);
397 37 : if (padfZ)
398 : {
399 18 : const double z = padfZ[i] + (padfZ[i + 1] - padfZ[i]) *
400 18 : (alpha - alpha0) /
401 18 : (alpha1 - alpha0);
402 18 : adfZ.emplace_back(z);
403 : }
404 37 : if (padfM)
405 : {
406 15 : adfM.emplace_back(padfM[i]);
407 : }
408 : }
409 : }
410 9 : aoRawPoint.emplace_back(x1, y1);
411 9 : if (padfZ)
412 4 : adfZ.emplace_back(padfZ[i + 1]);
413 9 : if (padfM)
414 3 : adfM.emplace_back(padfM[i + 1]);
415 :
416 9 : if (dfSegmentLength1 > dfMaxLength ||
417 : dfSegmentLength2 > dfMaxLength)
418 : {
419 9 : const double dfVal =
420 9 : 1 + 2 * std::floor(dfSegmentLength2 / dfMaxLength / 2.0);
421 18 : if (dfVal < 0.0 || std::isnan(dfVal) ||
422 9 : dfVal >= kMax - static_cast<int>(aoRawPoint.size()))
423 : {
424 0 : CPLError(CE_Failure, CPLE_AppDefined,
425 : "segmentize nIntermediatePoints invalid 2: %lf",
426 : dfVal);
427 0 : return false;
428 : }
429 9 : int nIntermediatePoints = static_cast<int>(dfVal);
430 9 : const double dfStep =
431 9 : (alpha2 - alpha1) / (nIntermediatePoints + 1);
432 46 : for (int j = 1; j <= nIntermediatePoints; ++j)
433 : {
434 37 : const double alpha = alpha1 + dfStep * j;
435 37 : const double x = cx + R * cos(alpha);
436 37 : const double y = cy + R * sin(alpha);
437 37 : aoRawPoint.emplace_back(x, y);
438 37 : if (padfZ)
439 : {
440 18 : const double z =
441 18 : padfZ[i + 1] + (padfZ[i + 2] - padfZ[i + 1]) *
442 18 : (alpha - alpha1) /
443 18 : (alpha2 - alpha1);
444 18 : adfZ.emplace_back(z);
445 : }
446 37 : if (padfM)
447 : {
448 15 : adfM.emplace_back(padfM[i + 1]);
449 : }
450 : }
451 : }
452 : }
453 : else
454 : {
455 : // It is a straight line.
456 2 : const double dfSegmentLength1 = dist(x0, y0, x1, y1);
457 2 : const double dfSegmentLength2 = dist(x1, y1, x2, y2);
458 2 : if (dfSegmentLength1 > dfMaxLength ||
459 : dfSegmentLength2 > dfMaxLength)
460 : {
461 2 : const double dfVal =
462 2 : 1 + 2 * std::ceil(dfSegmentLength1 / dfMaxLength / 2.0);
463 4 : if (dfVal < 0.0 || std::isnan(dfVal) ||
464 2 : dfVal >= kMax - static_cast<int>(aoRawPoint.size()))
465 : {
466 0 : CPLError(CE_Failure, CPLE_AppDefined,
467 : "segmentize nIntermediatePoints invalid 2: %lf",
468 : dfVal);
469 0 : return false;
470 : }
471 2 : int nIntermediatePoints = static_cast<int>(dfVal);
472 12 : for (int j = 1; j <= nIntermediatePoints; ++j)
473 : {
474 : aoRawPoint.emplace_back(
475 0 : x0 + j * (x1 - x0) / (nIntermediatePoints + 1),
476 10 : y0 + j * (y1 - y0) / (nIntermediatePoints + 1));
477 10 : if (padfZ)
478 5 : adfZ.emplace_back(padfZ[i] +
479 5 : j * (padfZ[i + 1] - padfZ[i]) /
480 5 : (nIntermediatePoints + 1));
481 10 : if (padfM)
482 0 : adfM.emplace_back(padfM[i]);
483 : }
484 : }
485 :
486 2 : aoRawPoint.emplace_back(x1, y1);
487 2 : if (padfZ)
488 1 : adfZ.emplace_back(padfZ[i + 1]);
489 2 : if (padfM)
490 0 : adfM.emplace_back(padfM[i + 1]);
491 :
492 2 : if (dfSegmentLength1 > dfMaxLength ||
493 : dfSegmentLength2 > dfMaxLength)
494 : {
495 2 : const double dfVal =
496 2 : 1 + 2 * std::ceil(dfSegmentLength2 / dfMaxLength / 2.0);
497 4 : if (dfVal < 0.0 || std::isnan(dfVal) ||
498 2 : dfVal >= kMax - static_cast<int>(aoRawPoint.size()))
499 : {
500 0 : CPLError(CE_Failure, CPLE_AppDefined,
501 : "segmentize nIntermediatePoints invalid 3: %lf",
502 : dfVal);
503 0 : return false;
504 : }
505 2 : const int nIntermediatePoints = static_cast<int>(dfVal);
506 :
507 12 : for (int j = 1; j <= nIntermediatePoints; ++j)
508 : {
509 : aoRawPoint.emplace_back(
510 0 : x1 + j * (x2 - x1) / (nIntermediatePoints + 1),
511 10 : y1 + j * (y2 - y1) / (nIntermediatePoints + 1));
512 10 : if (padfZ)
513 5 : adfZ.emplace_back(padfZ[i + 1] +
514 5 : j * (padfZ[i + 2] - padfZ[i + 1]) /
515 5 : (nIntermediatePoints + 1));
516 10 : if (padfM)
517 0 : adfM.emplace_back(padfM[i + 1]);
518 : }
519 : }
520 : }
521 : }
522 9 : aoRawPoint.push_back(paoPoints[nPointCount - 1]);
523 9 : if (padfZ)
524 4 : adfZ.push_back(padfZ[nPointCount - 1]);
525 9 : if (padfM)
526 2 : adfM.push_back(padfM[nPointCount - 1]);
527 :
528 9 : CPLAssert(aoRawPoint.empty() ||
529 : (aoRawPoint.size() >= 3 && (aoRawPoint.size() % 2) == 1));
530 9 : if (padfZ)
531 : {
532 4 : CPLAssert(adfZ.size() == aoRawPoint.size());
533 : }
534 9 : if (padfM)
535 : {
536 2 : CPLAssert(adfM.size() == aoRawPoint.size());
537 : }
538 :
539 : // Is there actually something to modify?
540 9 : if (nPointCount < static_cast<int>(aoRawPoint.size()))
541 : {
542 9 : nPointCount = static_cast<int>(aoRawPoint.size());
543 9 : paoPoints = static_cast<OGRRawPoint *>(
544 9 : CPLRealloc(paoPoints, sizeof(OGRRawPoint) * nPointCount));
545 9 : memcpy(paoPoints, &aoRawPoint[0], sizeof(OGRRawPoint) * nPointCount);
546 9 : if (padfZ)
547 : {
548 4 : padfZ = static_cast<double *>(
549 4 : CPLRealloc(padfZ, sizeof(double) * aoRawPoint.size()));
550 4 : memcpy(padfZ, &adfZ[0], sizeof(double) * nPointCount);
551 : }
552 9 : if (padfM)
553 : {
554 2 : padfM = static_cast<double *>(
555 2 : CPLRealloc(padfM, sizeof(double) * aoRawPoint.size()));
556 2 : memcpy(padfM, &adfM[0], sizeof(double) * nPointCount);
557 : }
558 : }
559 9 : return bRet;
560 : }
561 :
562 : /************************************************************************/
563 : /* Value() */
564 : /* */
565 : /* Get an interpolated point at some distance along the curve. */
566 : /************************************************************************/
567 :
568 14 : void OGRCircularString::Value(double dfDistance, OGRPoint *poPoint) const
569 :
570 : {
571 14 : if (dfDistance < 0)
572 : {
573 1 : StartPoint(poPoint);
574 1 : return;
575 : }
576 :
577 13 : double dfLength = 0;
578 :
579 19 : for (int i = 0; i < nPointCount - 2; i += 2)
580 : {
581 18 : const double x0 = paoPoints[i].x;
582 18 : const double y0 = paoPoints[i].y;
583 18 : const double x1 = paoPoints[i + 1].x;
584 18 : const double y1 = paoPoints[i + 1].y;
585 18 : const double x2 = paoPoints[i + 2].x;
586 18 : const double y2 = paoPoints[i + 2].y;
587 18 : double R = 0.0;
588 18 : double cx = 0.0;
589 18 : double cy = 0.0;
590 18 : double alpha0 = 0.0;
591 18 : double alpha1 = 0.0;
592 18 : double alpha2 = 0.0;
593 :
594 : // We have strong constraints on the number of intermediate points
595 : // we can add.
596 :
597 18 : if (OGRGeometryFactory::GetCurveParameters(
598 18 : x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
599 : {
600 : // It is an arc circle.
601 9 : const double dfSegLength = fabs(alpha2 - alpha0) * R;
602 9 : if (dfSegLength > 0)
603 : {
604 9 : if ((dfLength <= dfDistance) &&
605 9 : ((dfLength + dfSegLength) >= dfDistance))
606 : {
607 8 : const double dfRatio =
608 8 : (dfDistance - dfLength) / dfSegLength;
609 :
610 8 : const double alpha =
611 8 : alpha0 * (1 - dfRatio) + alpha2 * dfRatio;
612 8 : const double x = cx + R * cos(alpha);
613 8 : const double y = cy + R * sin(alpha);
614 :
615 8 : poPoint->setX(x);
616 8 : poPoint->setY(y);
617 :
618 8 : if (getCoordinateDimension() == 3)
619 2 : poPoint->setZ(padfZ[i] * (1 - dfRatio) +
620 2 : padfZ[i + 2] * dfRatio);
621 :
622 12 : return;
623 : }
624 :
625 1 : dfLength += dfSegLength;
626 : }
627 : }
628 : else
629 : {
630 : // It is a straight line.
631 9 : const double dfSegLength = dist(x0, y0, x2, y2);
632 9 : if (dfSegLength > 0)
633 : {
634 9 : if ((dfLength <= dfDistance) &&
635 9 : ((dfLength + dfSegLength) >= dfDistance))
636 : {
637 4 : const double dfRatio =
638 4 : (dfDistance - dfLength) / dfSegLength;
639 :
640 4 : poPoint->setX(paoPoints[i].x * (1 - dfRatio) +
641 4 : paoPoints[i + 2].x * dfRatio);
642 4 : poPoint->setY(paoPoints[i].y * (1 - dfRatio) +
643 4 : paoPoints[i + 2].y * dfRatio);
644 :
645 4 : if (getCoordinateDimension() == 3)
646 2 : poPoint->setZ(padfZ[i] * (1 - dfRatio) +
647 2 : padfZ[i + 2] * dfRatio);
648 :
649 4 : return;
650 : }
651 :
652 5 : dfLength += dfSegLength;
653 : }
654 : }
655 : }
656 :
657 1 : EndPoint(poPoint);
658 : }
659 :
660 : /************************************************************************/
661 : /* CurveToLine() */
662 : /************************************************************************/
663 :
664 : OGRLineString *
665 4077 : OGRCircularString::CurveToLine(double dfMaxAngleStepSizeDegrees,
666 : const char *const *papszOptions) const
667 : {
668 4077 : OGRLineString *poLine = new OGRLineString();
669 4077 : poLine->assignSpatialReference(getSpatialReference());
670 :
671 4077 : const bool bHasZ = getCoordinateDimension() == 3;
672 8371 : for (int i = 0; i < nPointCount - 2; i += 2)
673 : {
674 12066 : OGRLineString *poArc = OGRGeometryFactory::curveToLineString(
675 4294 : paoPoints[i].x, paoPoints[i].y, padfZ ? padfZ[i] : 0.0,
676 4294 : paoPoints[i + 1].x, paoPoints[i + 1].y, padfZ ? padfZ[i + 1] : 0.0,
677 8316 : paoPoints[i + 2].x, paoPoints[i + 2].y, padfZ ? padfZ[i + 2] : 0.0,
678 : bHasZ, dfMaxAngleStepSizeDegrees, papszOptions);
679 4294 : poLine->addSubLineString(poArc, (i == 0) ? 0 : 1);
680 4294 : delete poArc;
681 : }
682 4077 : return poLine;
683 : }
684 :
685 : /************************************************************************/
686 : /* IsValidFast() */
687 : /************************************************************************/
688 :
689 4206 : bool OGRCircularString::IsValidFast(std::string *posReason) const
690 :
691 : {
692 4206 : if (nPointCount == 1 || nPointCount == 2 ||
693 4205 : (nPointCount >= 3 && (nPointCount % 2) == 0))
694 : {
695 1 : if (posReason)
696 : {
697 : *posReason =
698 : CPLSPrintf("Invalid number of points in circular string : %d",
699 0 : nPointCount);
700 : }
701 : else
702 : {
703 1 : CPLError(CE_Failure, CPLE_NotSupported,
704 : "Invalid number of points in circular string : %d",
705 1 : nPointCount);
706 : }
707 1 : return FALSE;
708 : }
709 4205 : return TRUE;
710 : }
711 :
712 : /************************************************************************/
713 : /* IsValid() */
714 : /************************************************************************/
715 :
716 2 : bool OGRCircularString::IsValid(std::string *posReason) const
717 :
718 : {
719 2 : return IsValidFast(posReason) && OGRGeometry::IsValid(posReason);
720 : }
721 :
722 : /************************************************************************/
723 : /* hasCurveGeometry() */
724 : /************************************************************************/
725 :
726 3974 : bool OGRCircularString::hasCurveGeometry(int /* bLookForNonLinear */) const
727 : {
728 3974 : return TRUE;
729 : }
730 :
731 : /************************************************************************/
732 : /* getLinearGeometry() */
733 : /************************************************************************/
734 :
735 : OGRGeometry *
736 3310 : OGRCircularString::getLinearGeometry(double dfMaxAngleStepSizeDegrees,
737 : const char *const *papszOptions) const
738 : {
739 3310 : return CurveToLine(dfMaxAngleStepSizeDegrees, papszOptions);
740 : }
741 :
742 : //! @cond Doxygen_Suppress
743 : /************************************************************************/
744 : /* GetCasterToLineString() */
745 : /************************************************************************/
746 :
747 0 : static OGRLineString *CasterToLineString(OGRCurve *poGeom)
748 : {
749 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible",
750 0 : poGeom->getGeometryName());
751 0 : delete poGeom;
752 0 : return nullptr;
753 : }
754 :
755 0 : OGRCurveCasterToLineString OGRCircularString::GetCasterToLineString() const
756 : {
757 0 : return ::CasterToLineString;
758 : }
759 :
760 : /************************************************************************/
761 : /* GetCasterToLinearRing() */
762 : /************************************************************************/
763 :
764 0 : static OGRLinearRing *CasterToLinearRing(OGRCurve *poGeom)
765 : {
766 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible",
767 0 : poGeom->getGeometryName());
768 0 : delete poGeom;
769 0 : return nullptr;
770 : }
771 :
772 0 : OGRCurveCasterToLinearRing OGRCircularString::GetCasterToLinearRing() const
773 : {
774 0 : return ::CasterToLinearRing;
775 : }
776 :
777 : //! @endcond
778 :
779 : /************************************************************************/
780 : /* IsFullCircle() */
781 : /************************************************************************/
782 :
783 12 : int OGRCircularString::IsFullCircle(double &cx, double &cy,
784 : double &square_R) const
785 : {
786 12 : if (getNumPoints() == 3 && get_IsClosed())
787 : {
788 7 : const double x0 = getX(0);
789 7 : const double y0 = getY(0);
790 7 : const double x1 = getX(1);
791 7 : const double y1 = getY(1);
792 7 : cx = (x0 + x1) / 2;
793 7 : cy = (y0 + y1) / 2;
794 7 : square_R = (x1 - cx) * (x1 - cx) + (y1 - cy) * (y1 - cy);
795 7 : return TRUE;
796 : }
797 : // Full circle defined by 2 arcs?
798 5 : else if (getNumPoints() == 5 && get_IsClosed())
799 : {
800 4 : double R_1 = 0.0;
801 4 : double cx_1 = 0.0;
802 4 : double cy_1 = 0.0;
803 4 : double alpha0_1 = 0.0;
804 4 : double alpha1_1 = 0.0;
805 4 : double alpha2_1 = 0.0;
806 4 : double R_2 = 0.0;
807 4 : double cx_2 = 0.0;
808 4 : double cy_2 = 0.0;
809 4 : double alpha0_2 = 0.0;
810 4 : double alpha1_2 = 0.0;
811 4 : double alpha2_2 = 0.0;
812 4 : if (OGRGeometryFactory::GetCurveParameters(
813 : getX(0), getY(0), getX(1), getY(1), getX(2), getY(2), R_1, cx_1,
814 4 : cy_1, alpha0_1, alpha1_1, alpha2_1) &&
815 4 : OGRGeometryFactory::GetCurveParameters(
816 : getX(2), getY(2), getX(3), getY(3), getX(4), getY(4), R_2, cx_2,
817 4 : cy_2, alpha0_2, alpha1_2, alpha2_2) &&
818 4 : fabs(R_1 - R_2) < 1e-10 && fabs(cx_1 - cx_2) < 1e-10 &&
819 12 : fabs(cy_1 - cy_2) < 1e-10 &&
820 4 : (alpha2_1 - alpha0_1) * (alpha2_2 - alpha0_2) > 0)
821 : {
822 3 : cx = cx_1;
823 3 : cy = cy_1;
824 3 : square_R = R_1 * R_1;
825 3 : return TRUE;
826 : }
827 : }
828 2 : return FALSE;
829 : }
830 :
831 : /************************************************************************/
832 : /* get_AreaOfCurveSegments() */
833 : /************************************************************************/
834 :
835 : //! @cond Doxygen_Suppress
836 25 : double OGRCircularString::get_AreaOfCurveSegments() const
837 : {
838 25 : double dfArea = 0.0;
839 53 : for (int i = 0; i < getNumPoints() - 2; i += 2)
840 : {
841 28 : const double x0 = getX(i);
842 28 : const double y0 = getY(i);
843 28 : const double x1 = getX(i + 1);
844 28 : const double y1 = getY(i + 1);
845 28 : const double x2 = getX(i + 2);
846 28 : const double y2 = getY(i + 2);
847 28 : double R = 0.0;
848 28 : double cx = 0.0;
849 28 : double cy = 0.0;
850 28 : double alpha0 = 0.0;
851 28 : double alpha1 = 0.0;
852 28 : double alpha2 = 0.0;
853 28 : if (OGRGeometryFactory::GetCurveParameters(
854 28 : x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
855 : {
856 : // Should be <= PI in absolute value.
857 26 : const double delta_alpha01 = alpha1 - alpha0;
858 26 : const double delta_alpha12 = alpha2 - alpha1; // Same.
859 : // http://en.wikipedia.org/wiki/Circular_segment
860 26 : dfArea += 0.5 * R * R *
861 26 : fabs(delta_alpha01 - sin(delta_alpha01) + delta_alpha12 -
862 26 : sin(delta_alpha12));
863 : }
864 : }
865 25 : return dfArea;
866 : }
867 :
868 : //! @endcond
869 :
870 : /************************************************************************/
871 : /* get_Area() */
872 : /************************************************************************/
873 :
874 5 : double OGRCircularString::get_Area() const
875 : {
876 5 : if (IsEmpty() || !get_IsClosed())
877 1 : return 0;
878 :
879 4 : double cx = 0.0;
880 4 : double cy = 0.0;
881 4 : double square_R = 0.0;
882 :
883 4 : if (IsFullCircle(cx, cy, square_R))
884 : {
885 3 : return M_PI * square_R;
886 : }
887 :
888 : // Optimization for convex rings.
889 1 : if (IsConvex())
890 : {
891 : // Compute area of shape without the circular segments.
892 1 : double dfArea = get_LinearArea();
893 :
894 : // Add the area of the spherical segments.
895 1 : dfArea += get_AreaOfCurveSegments();
896 :
897 1 : return dfArea;
898 : }
899 :
900 0 : OGRLineString *poLS = CurveToLine();
901 0 : const double dfArea = poLS->get_Area();
902 0 : delete poLS;
903 :
904 0 : return dfArea;
905 : }
906 :
907 : /************************************************************************/
908 : /* get_GeodesicArea() */
909 : /************************************************************************/
910 :
911 4 : double OGRCircularString::get_GeodesicArea(
912 : const OGRSpatialReference *poSRSOverride) const
913 : {
914 4 : if (IsEmpty())
915 1 : return 0;
916 :
917 3 : if (!get_IsClosed())
918 : {
919 2 : CPLError(CE_Failure, CPLE_AppDefined, "Non-closed geometry");
920 2 : return -1;
921 : }
922 :
923 1 : if (!poSRSOverride)
924 1 : poSRSOverride = getSpatialReference();
925 :
926 2 : auto poLS = std::unique_ptr<OGRLineString>(CurveToLine());
927 1 : return poLS->get_GeodesicArea(poSRSOverride);
928 : }
929 :
930 : /************************************************************************/
931 : /* get_GeodesicLength() */
932 : /************************************************************************/
933 :
934 2 : double OGRCircularString::get_GeodesicLength(
935 : const OGRSpatialReference *poSRSOverride) const
936 : {
937 2 : if (IsEmpty())
938 1 : return 0;
939 :
940 1 : if (!poSRSOverride)
941 1 : poSRSOverride = getSpatialReference();
942 :
943 2 : auto poLS = std::unique_ptr<OGRLineString>(CurveToLine());
944 1 : return poLS->get_GeodesicLength(poSRSOverride);
945 : }
946 :
947 : //! @cond Doxygen_Suppress
948 :
949 : /************************************************************************/
950 : /* ContainsPoint() */
951 : /************************************************************************/
952 :
953 6 : int OGRCircularString::ContainsPoint(const OGRPoint *p) const
954 : {
955 6 : double cx = 0.0;
956 6 : double cy = 0.0;
957 6 : double square_R = 0.0;
958 6 : if (IsFullCircle(cx, cy, square_R))
959 : {
960 5 : const double square_dist = (p->getX() - cx) * (p->getX() - cx) +
961 5 : (p->getY() - cy) * (p->getY() - cy);
962 5 : return square_dist < square_R;
963 : }
964 1 : return -1;
965 : }
966 :
967 : /************************************************************************/
968 : /* IntersectsPoint() */
969 : /************************************************************************/
970 :
971 2 : int OGRCircularString::IntersectsPoint(const OGRPoint *p) const
972 : {
973 2 : double cx = 0.0;
974 2 : double cy = 0.0;
975 2 : double square_R = 0.0;
976 2 : if (IsFullCircle(cx, cy, square_R))
977 : {
978 2 : const double square_dist = (p->getX() - cx) * (p->getX() - cx) +
979 2 : (p->getY() - cy) * (p->getY() - cy);
980 2 : return square_dist <= square_R;
981 : }
982 0 : return -1;
983 : }
984 :
985 : //! @endcond
|
