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