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 10549 : OGRwkbGeometryType OGRCircularString::getGeometryType() const
66 :
67 : {
68 10549 : if ((flags & OGR_G_3D) && (flags & OGR_G_MEASURED))
69 683 : return wkbCircularStringZM;
70 9866 : else if (flags & OGR_G_MEASURED)
71 67 : return wkbCircularStringM;
72 9799 : else if (flags & OGR_G_3D)
73 539 : return wkbCircularStringZ;
74 : else
75 9260 : return wkbCircularString;
76 : }
77 :
78 : /************************************************************************/
79 : /* getGeometryName() */
80 : /************************************************************************/
81 :
82 4839 : const char *OGRCircularString::getGeometryName() const
83 :
84 : {
85 4839 : return "CIRCULARSTRING";
86 : }
87 :
88 : /************************************************************************/
89 : /* clone() */
90 : /************************************************************************/
91 :
92 503 : OGRCircularString *OGRCircularString::clone() const
93 :
94 : {
95 503 : 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 3550 : OGRErr OGRCircularString::importFromWkt(const char **ppszInput)
156 :
157 : {
158 3550 : const OGRErr eErr = OGRSimpleCurve::importFromWkt(ppszInput);
159 3550 : if (eErr == OGRERR_NONE)
160 : {
161 3550 : if (!IsValidFast())
162 : {
163 1 : empty();
164 1 : return OGRERR_CORRUPT_DATA;
165 : }
166 : }
167 3549 : 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 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 4058 : OGRCircularString::CurveToLine(double dfMaxAngleStepSizeDegrees,
666 : const char *const *papszOptions) const
667 : {
668 4058 : OGRLineString *poLine = new OGRLineString();
669 4058 : poLine->assignSpatialReference(getSpatialReference());
670 :
671 4058 : const bool bHasZ = getCoordinateDimension() == 3;
672 8333 : for (int i = 0; i < nPointCount - 2; i += 2)
673 : {
674 12009 : OGRLineString *poArc = OGRGeometryFactory::curveToLineString(
675 4275 : paoPoints[i].x, paoPoints[i].y, padfZ ? padfZ[i] : 0.0,
676 4275 : paoPoints[i + 1].x, paoPoints[i + 1].y, padfZ ? padfZ[i + 1] : 0.0,
677 8278 : paoPoints[i + 2].x, paoPoints[i + 2].y, padfZ ? padfZ[i + 2] : 0.0,
678 : bHasZ, dfMaxAngleStepSizeDegrees, papszOptions);
679 4275 : poLine->addSubLineString(poArc, (i == 0) ? 0 : 1);
680 4275 : delete poArc;
681 : }
682 4058 : return poLine;
683 : }
684 :
685 : /************************************************************************/
686 : /* IsValidFast() */
687 : /************************************************************************/
688 :
689 4192 : OGRBoolean OGRCircularString::IsValidFast() const
690 :
691 : {
692 4192 : if (nPointCount == 1 || nPointCount == 2 ||
693 4191 : (nPointCount >= 3 && (nPointCount % 2) == 0))
694 : {
695 1 : CPLError(CE_Failure, CPLE_NotSupported,
696 1 : "Bad number of points in circular string : %d", nPointCount);
697 1 : return FALSE;
698 : }
699 4191 : return TRUE;
700 : }
701 :
702 : /************************************************************************/
703 : /* IsValid() */
704 : /************************************************************************/
705 :
706 2 : OGRBoolean OGRCircularString::IsValid() const
707 :
708 : {
709 2 : return IsValidFast() && OGRGeometry::IsValid();
710 : }
711 :
712 : /************************************************************************/
713 : /* hasCurveGeometry() */
714 : /************************************************************************/
715 :
716 : OGRBoolean
717 3968 : OGRCircularString::hasCurveGeometry(int /* bLookForNonLinear */) const
718 : {
719 3968 : return TRUE;
720 : }
721 :
722 : /************************************************************************/
723 : /* getLinearGeometry() */
724 : /************************************************************************/
725 :
726 : OGRGeometry *
727 3312 : OGRCircularString::getLinearGeometry(double dfMaxAngleStepSizeDegrees,
728 : const char *const *papszOptions) const
729 : {
730 3312 : return CurveToLine(dfMaxAngleStepSizeDegrees, papszOptions);
731 : }
732 :
733 : //! @cond Doxygen_Suppress
734 : /************************************************************************/
735 : /* GetCasterToLineString() */
736 : /************************************************************************/
737 :
738 0 : static OGRLineString *CasterToLineString(OGRCurve *poGeom)
739 : {
740 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible",
741 0 : poGeom->getGeometryName());
742 0 : delete poGeom;
743 0 : return nullptr;
744 : }
745 :
746 0 : OGRCurveCasterToLineString OGRCircularString::GetCasterToLineString() const
747 : {
748 0 : return ::CasterToLineString;
749 : }
750 :
751 : /************************************************************************/
752 : /* GetCasterToLinearRing() */
753 : /************************************************************************/
754 :
755 0 : static OGRLinearRing *CasterToLinearRing(OGRCurve *poGeom)
756 : {
757 0 : CPLError(CE_Failure, CPLE_AppDefined, "%s found. Conversion impossible",
758 0 : poGeom->getGeometryName());
759 0 : delete poGeom;
760 0 : return nullptr;
761 : }
762 :
763 0 : OGRCurveCasterToLinearRing OGRCircularString::GetCasterToLinearRing() const
764 : {
765 0 : return ::CasterToLinearRing;
766 : }
767 :
768 : //! @endcond
769 :
770 : /************************************************************************/
771 : /* IsFullCircle() */
772 : /************************************************************************/
773 :
774 12 : int OGRCircularString::IsFullCircle(double &cx, double &cy,
775 : double &square_R) const
776 : {
777 12 : if (getNumPoints() == 3 && get_IsClosed())
778 : {
779 7 : const double x0 = getX(0);
780 7 : const double y0 = getY(0);
781 7 : const double x1 = getX(1);
782 7 : const double y1 = getY(1);
783 7 : cx = (x0 + x1) / 2;
784 7 : cy = (y0 + y1) / 2;
785 7 : square_R = (x1 - cx) * (x1 - cx) + (y1 - cy) * (y1 - cy);
786 7 : return TRUE;
787 : }
788 : // Full circle defined by 2 arcs?
789 5 : else if (getNumPoints() == 5 && get_IsClosed())
790 : {
791 4 : double R_1 = 0.0;
792 4 : double cx_1 = 0.0;
793 4 : double cy_1 = 0.0;
794 4 : double alpha0_1 = 0.0;
795 4 : double alpha1_1 = 0.0;
796 4 : double alpha2_1 = 0.0;
797 4 : double R_2 = 0.0;
798 4 : double cx_2 = 0.0;
799 4 : double cy_2 = 0.0;
800 4 : double alpha0_2 = 0.0;
801 4 : double alpha1_2 = 0.0;
802 4 : double alpha2_2 = 0.0;
803 4 : if (OGRGeometryFactory::GetCurveParameters(
804 : getX(0), getY(0), getX(1), getY(1), getX(2), getY(2), R_1, cx_1,
805 4 : cy_1, alpha0_1, alpha1_1, alpha2_1) &&
806 4 : OGRGeometryFactory::GetCurveParameters(
807 : getX(2), getY(2), getX(3), getY(3), getX(4), getY(4), R_2, cx_2,
808 4 : cy_2, alpha0_2, alpha1_2, alpha2_2) &&
809 4 : fabs(R_1 - R_2) < 1e-10 && fabs(cx_1 - cx_2) < 1e-10 &&
810 12 : fabs(cy_1 - cy_2) < 1e-10 &&
811 4 : (alpha2_1 - alpha0_1) * (alpha2_2 - alpha0_2) > 0)
812 : {
813 3 : cx = cx_1;
814 3 : cy = cy_1;
815 3 : square_R = R_1 * R_1;
816 3 : return TRUE;
817 : }
818 : }
819 2 : return FALSE;
820 : }
821 :
822 : /************************************************************************/
823 : /* get_AreaOfCurveSegments() */
824 : /************************************************************************/
825 :
826 : //! @cond Doxygen_Suppress
827 25 : double OGRCircularString::get_AreaOfCurveSegments() const
828 : {
829 25 : double dfArea = 0.0;
830 53 : for (int i = 0; i < getNumPoints() - 2; i += 2)
831 : {
832 28 : const double x0 = getX(i);
833 28 : const double y0 = getY(i);
834 28 : const double x1 = getX(i + 1);
835 28 : const double y1 = getY(i + 1);
836 28 : const double x2 = getX(i + 2);
837 28 : const double y2 = getY(i + 2);
838 28 : double R = 0.0;
839 28 : double cx = 0.0;
840 28 : double cy = 0.0;
841 28 : double alpha0 = 0.0;
842 28 : double alpha1 = 0.0;
843 28 : double alpha2 = 0.0;
844 28 : if (OGRGeometryFactory::GetCurveParameters(
845 28 : x0, y0, x1, y1, x2, y2, R, cx, cy, alpha0, alpha1, alpha2))
846 : {
847 : // Should be <= PI in absolute value.
848 26 : const double delta_alpha01 = alpha1 - alpha0;
849 26 : const double delta_alpha12 = alpha2 - alpha1; // Same.
850 : // http://en.wikipedia.org/wiki/Circular_segment
851 26 : dfArea += 0.5 * R * R *
852 26 : fabs(delta_alpha01 - sin(delta_alpha01) + delta_alpha12 -
853 26 : sin(delta_alpha12));
854 : }
855 : }
856 25 : return dfArea;
857 : }
858 :
859 : //! @endcond
860 :
861 : /************************************************************************/
862 : /* get_Area() */
863 : /************************************************************************/
864 :
865 5 : double OGRCircularString::get_Area() const
866 : {
867 5 : if (IsEmpty() || !get_IsClosed())
868 1 : return 0;
869 :
870 4 : double cx = 0.0;
871 4 : double cy = 0.0;
872 4 : double square_R = 0.0;
873 :
874 4 : if (IsFullCircle(cx, cy, square_R))
875 : {
876 3 : return M_PI * square_R;
877 : }
878 :
879 : // Optimization for convex rings.
880 1 : if (IsConvex())
881 : {
882 : // Compute area of shape without the circular segments.
883 1 : double dfArea = get_LinearArea();
884 :
885 : // Add the area of the spherical segments.
886 1 : dfArea += get_AreaOfCurveSegments();
887 :
888 1 : return dfArea;
889 : }
890 :
891 0 : OGRLineString *poLS = CurveToLine();
892 0 : const double dfArea = poLS->get_Area();
893 0 : delete poLS;
894 :
895 0 : return dfArea;
896 : }
897 :
898 : /************************************************************************/
899 : /* get_GeodesicArea() */
900 : /************************************************************************/
901 :
902 4 : double OGRCircularString::get_GeodesicArea(
903 : const OGRSpatialReference *poSRSOverride) const
904 : {
905 4 : if (IsEmpty())
906 1 : return 0;
907 :
908 3 : if (!get_IsClosed())
909 : {
910 2 : CPLError(CE_Failure, CPLE_AppDefined, "Non-closed geometry");
911 2 : return -1;
912 : }
913 :
914 1 : if (!poSRSOverride)
915 1 : poSRSOverride = getSpatialReference();
916 :
917 2 : auto poLS = std::unique_ptr<OGRLineString>(CurveToLine());
918 1 : return poLS->get_GeodesicArea(poSRSOverride);
919 : }
920 :
921 : /************************************************************************/
922 : /* get_GeodesicLength() */
923 : /************************************************************************/
924 :
925 2 : double OGRCircularString::get_GeodesicLength(
926 : const OGRSpatialReference *poSRSOverride) const
927 : {
928 2 : if (IsEmpty())
929 1 : return 0;
930 :
931 1 : if (!poSRSOverride)
932 1 : poSRSOverride = getSpatialReference();
933 :
934 2 : auto poLS = std::unique_ptr<OGRLineString>(CurveToLine());
935 1 : return poLS->get_GeodesicLength(poSRSOverride);
936 : }
937 :
938 : //! @cond Doxygen_Suppress
939 :
940 : /************************************************************************/
941 : /* ContainsPoint() */
942 : /************************************************************************/
943 :
944 6 : int OGRCircularString::ContainsPoint(const OGRPoint *p) const
945 : {
946 6 : double cx = 0.0;
947 6 : double cy = 0.0;
948 6 : double square_R = 0.0;
949 6 : if (IsFullCircle(cx, cy, square_R))
950 : {
951 5 : const double square_dist = (p->getX() - cx) * (p->getX() - cx) +
952 5 : (p->getY() - cy) * (p->getY() - cy);
953 5 : return square_dist < square_R;
954 : }
955 1 : return -1;
956 : }
957 :
958 : /************************************************************************/
959 : /* IntersectsPoint() */
960 : /************************************************************************/
961 :
962 2 : int OGRCircularString::IntersectsPoint(const OGRPoint *p) const
963 : {
964 2 : double cx = 0.0;
965 2 : double cy = 0.0;
966 2 : double square_R = 0.0;
967 2 : if (IsFullCircle(cx, cy, square_R))
968 : {
969 2 : const double square_dist = (p->getX() - cx) * (p->getX() - cx) +
970 2 : (p->getY() - cy) * (p->getY() - cy);
971 2 : return square_dist <= square_R;
972 : }
973 0 : return -1;
974 : }
975 :
976 : //! @endcond
|
