LCOV - code coverage report
Current view: top level - ogr - ogrcircularstring.cpp (source / functions) Hit Total Coverage
Test: gdal_filtered.info Lines: 406 452 89.8 %
Date: 2025-10-21 22:35:35 Functions: 27 31 87.1 %

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

Generated by: LCOV version 1.14