/src/gdal/ogr/ogrlinearring.cpp

Source (jump to first uncovered line)
/******************************************************************************
 *
 * Project:  OpenGIS Simple Features Reference Implementation
 * Purpose:  The OGRLinearRing geometry class.
 * Author:   Frank Warmerdam, warmerdam@pobox.com
 *
 ******************************************************************************
 * Copyright (c) 1999, Frank Warmerdam
 * Copyright (c) 2008-2014, Even Rouault <even dot rouault at spatialys.com>
 *
 * SPDX-License-Identifier: MIT
 ****************************************************************************/

#include "cpl_port.h"
#include "ogr_geometry.h"

#include <climits>
#include <cmath>
#include <cstring>
#include <limits>

#include "cpl_error.h"
#include "ogr_core.h"
#include "ogr_geometry.h"
#include "ogr_p.h"

/************************************************************************/
/*                  OGRLinearRing( const OGRLinearRing& )               */
/************************************************************************/

/**
 * \brief Copy constructor.
 *
 * Note: before GDAL 2.1, only the default implementation of the constructor
 * existed, which could be unsafe to use.
 *
 * @since GDAL 2.1
 */

OGRLinearRing::OGRLinearRing(const OGRLinearRing &) = default;

/************************************************************************/
/*                           OGRLinearRing()                            */
/************************************************************************/

/** Constructor
 * @param poSrcRing source ring.
 */
OGRLinearRing::OGRLinearRing(const OGRLinearRing *poSrcRing)

{
    if (poSrcRing == nullptr)
    {
        CPLDebug("OGR",
                 "OGRLinearRing::OGRLinearRing(OGRLinearRing*poSrcRing) - "
                 "passed in ring is NULL!");
        return;
    }

    setNumPoints(poSrcRing->getNumPoints(), FALSE);

    memcpy(paoPoints, poSrcRing->paoPoints,
           sizeof(OGRRawPoint) * getNumPoints());

    if (poSrcRing->padfZ)
    {
        Make3D();
        memcpy(padfZ, poSrcRing->padfZ, sizeof(double) * getNumPoints());
    }
}

/************************************************************************/
/*                    operator=( const OGRLinearRing& )                 */
/************************************************************************/

/**
 * \brief Assignment operator.
 *
 * Note: before GDAL 2.1, only the default implementation of the operator
 * existed, which could be unsafe to use.
 *
 * @since GDAL 2.1
 */

OGRLinearRing &OGRLinearRing::operator=(const OGRLinearRing &other)
{
    if (this != &other)
    {
        OGRLineString::operator=(other);
    }
    return *this;
}

/************************************************************************/
/*                          getGeometryName()                           */
/************************************************************************/

const char *OGRLinearRing::getGeometryName() const

{
    return "LINEARRING";
}

/************************************************************************/
/*                              WkbSize()                               */
/*                                                                      */
/*      Disable this method.                                            */
/************************************************************************/

size_t OGRLinearRing::WkbSize() const

{
    return 0;
}

/************************************************************************/
/*                           importFromWkb()                            */
/*                                                                      */
/*      Disable method for this class.                                  */
/************************************************************************/

OGRErr OGRLinearRing::importFromWkb(const unsigned char * /*pabyData*/,
                                    size_t /*nSize*/,
                                    OGRwkbVariant /*eWkbVariant*/,
                                    size_t & /* nBytesConsumedOut */)

{
    return OGRERR_UNSUPPORTED_OPERATION;
}

/************************************************************************/
/*                            exportToWkb()                             */
/*                                                                      */
/*      Disable method for this class.                                  */
/************************************************************************/

OGRErr OGRLinearRing::exportToWkb(CPL_UNUSED unsigned char *pabyData,
                                  CPL_UNUSED const OGRwkbExportOptions *) const

{
    return OGRERR_UNSUPPORTED_OPERATION;
}

/************************************************************************/
/*                           _importFromWkb()                           */
/*                                                                      */
/*      Helper method for OGRPolygon.  NOT A NORMAL importFromWkb()     */
/*      method.                                                         */
/************************************************************************/

//! @cond Doxygen_Suppress
OGRErr OGRLinearRing::_importFromWkb(OGRwkbByteOrder eByteOrder, int _flags,
                                     const unsigned char *pabyData,
                                     size_t nBytesAvailable,
                                     size_t &nBytesConsumedOut)

{
    nBytesConsumedOut = 0;
    if (nBytesAvailable < 4 && nBytesAvailable != static_cast<size_t>(-1))
        return OGRERR_NOT_ENOUGH_DATA;

    /* -------------------------------------------------------------------- */
    /*      Get the vertex count.                                           */
    /* -------------------------------------------------------------------- */
    int nNewNumPoints = 0;

    memcpy(&nNewNumPoints, pabyData, 4);

    if (OGR_SWAP(eByteOrder))
        nNewNumPoints = CPL_SWAP32(nNewNumPoints);

    // Check if the wkb stream buffer is big enough to store
    // fetched number of points.
    // 16, 24, or 32 - size of point structure.
    size_t nPointSize = 0;
    if ((_flags & OGR_G_3D) && (_flags & OGR_G_MEASURED))
        nPointSize = 32;
    else if ((_flags & OGR_G_3D) || (_flags & OGR_G_MEASURED))
        nPointSize = 24;
    else
        nPointSize = 16;

    if (nNewNumPoints < 0 ||
        static_cast<size_t>(nNewNumPoints) >
            std::numeric_limits<size_t>::max() / nPointSize)
    {
        return OGRERR_CORRUPT_DATA;
    }
    const size_t nBufferMinSize = nPointSize * nNewNumPoints;
    if (nBytesAvailable != static_cast<size_t>(-1) &&
        nBufferMinSize > nBytesAvailable - 4)
    {
        CPLError(CE_Failure, CPLE_AppDefined,
                 "Length of input WKB is too small");
        return OGRERR_NOT_ENOUGH_DATA;
    }

    // (Re)Allocation of paoPoints buffer.
    setNumPoints(nNewNumPoints, FALSE);

    if (_flags & OGR_G_3D)
        Make3D();
    else
        Make2D();

    if (_flags & OGR_G_MEASURED)
        AddM();
    else
        RemoveM();

    nBytesConsumedOut = 4 + nPointCount * nPointSize;

    /* -------------------------------------------------------------------- */
    /*      Get the vertices                                                */
    /* -------------------------------------------------------------------- */
    if ((flags & OGR_G_3D) && (flags & OGR_G_MEASURED))
    {
        for (size_t i = 0; i < static_cast<size_t>(nPointCount); i++)
        {
            memcpy(&(paoPoints[i].x), pabyData + 4 + 32 * i, 8);
            memcpy(&(paoPoints[i].y), pabyData + 4 + 32 * i + 8, 8);
            memcpy(padfZ + i, pabyData + 4 + 32 * i + 16, 8);
            memcpy(padfM + i, pabyData + 4 + 32 * i + 24, 8);
        }
    }
    else if (flags & OGR_G_MEASURED)
    {
        for (size_t i = 0; i < static_cast<size_t>(nPointCount); i++)
        {
            memcpy(&(paoPoints[i].x), pabyData + 4 + 24 * i, 8);
            memcpy(&(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8);
            memcpy(padfM + i, pabyData + 4 + 24 * i + 16, 8);
        }
    }
    else if (flags & OGR_G_3D)
    {
        for (size_t i = 0; i < static_cast<size_t>(nPointCount); i++)
        {
            memcpy(&(paoPoints[i].x), pabyData + 4 + 24 * i, 8);
            memcpy(&(paoPoints[i].y), pabyData + 4 + 24 * i + 8, 8);
            memcpy(padfZ + i, pabyData + 4 + 24 * i + 16, 8);
        }
    }
    else if (nPointCount != 0)
    {
        memcpy(paoPoints, pabyData + 4, 16 * static_cast<size_t>(nPointCount));
    }

    /* -------------------------------------------------------------------- */
    /*      Byte swap if needed.                                            */
    /* -------------------------------------------------------------------- */
    if (OGR_SWAP(eByteOrder))
    {
        for (size_t i = 0; i < static_cast<size_t>(nPointCount); i++)
        {
            CPL_SWAPDOUBLE(&(paoPoints[i].x));
            CPL_SWAPDOUBLE(&(paoPoints[i].y));

            if (flags & OGR_G_3D)
            {
                CPL_SWAPDOUBLE(padfZ + i);
            }
            if (flags & OGR_G_MEASURED)
            {
                CPL_SWAPDOUBLE(padfM + i);
            }
        }
    }

    return OGRERR_NONE;
}

/************************************************************************/
/*                            _exportToWkb()                            */
/*                                                                      */
/*      Helper method for OGRPolygon.  THIS IS NOT THE NORMAL           */
/*      exportToWkb() METHOD.                                           */
/************************************************************************/

OGRErr OGRLinearRing::_exportToWkb(int _flags, unsigned char *pabyData,
                                   const OGRwkbExportOptions *psOptions) const

{

    /* -------------------------------------------------------------------- */
    /*      Copy in the raw data.                                           */
    /* -------------------------------------------------------------------- */
    memcpy(pabyData, &nPointCount, 4);

    /* -------------------------------------------------------------------- */
    /*      Copy in the raw data.                                           */
    /* -------------------------------------------------------------------- */
    size_t nWords = 0;
    if ((_flags & OGR_G_3D) && (_flags & OGR_G_MEASURED))
    {
        nWords = 4 * static_cast<size_t>(nPointCount);
        for (size_t i = 0; i < static_cast<size_t>(nPointCount); i++)
        {
            memcpy(pabyData + 4 + i * 32, &(paoPoints[i].x), 8);
            memcpy(pabyData + 4 + i * 32 + 8, &(paoPoints[i].y), 8);
            if (padfZ == nullptr)
                memset(pabyData + 4 + i * 32 + 16, 0, 8);
            else
                memcpy(pabyData + 4 + i * 32 + 16, padfZ + i, 8);
            if (padfM == nullptr)
                memset(pabyData + 4 + i * 32 + 24, 0, 8);
            else
                memcpy(pabyData + 4 + i * 32 + 24, padfM + i, 8);
        }
        OGRRoundCoordinatesIEEE754XYValues<32>(
            psOptions->sPrecision.nXYBitPrecision, pabyData + 4, nPointCount);
        OGRRoundCoordinatesIEEE754<32>(psOptions->sPrecision.nZBitPrecision,
                                       pabyData + 4 + 2 * sizeof(uint64_t),
                                       nPointCount);
        OGRRoundCoordinatesIEEE754<32>(psOptions->sPrecision.nMBitPrecision,
                                       pabyData + 4 + 3 * sizeof(uint64_t),
                                       nPointCount);
    }
    else if (_flags & OGR_G_MEASURED)
    {
        nWords = 3 * static_cast<size_t>(nPointCount);
        for (size_t i = 0; i < static_cast<size_t>(nPointCount); i++)
        {
            memcpy(pabyData + 4 + i * 24, &(paoPoints[i].x), 8);
            memcpy(pabyData + 4 + i * 24 + 8, &(paoPoints[i].y), 8);
            if (padfM == nullptr)
                memset(pabyData + 4 + i * 24 + 16, 0, 8);
            else
                memcpy(pabyData + 4 + i * 24 + 16, padfM + i, 8);
        }
        OGRRoundCoordinatesIEEE754XYValues<24>(
            psOptions->sPrecision.nXYBitPrecision, pabyData + 4, nPointCount);
        OGRRoundCoordinatesIEEE754<24>(psOptions->sPrecision.nMBitPrecision,
                                       pabyData + 4 + 2 * sizeof(uint64_t),
                                       nPointCount);
    }
    else if (_flags & OGR_G_3D)
    {
        nWords = 3 * static_cast<size_t>(nPointCount);
        for (size_t i = 0; i < static_cast<size_t>(nPointCount); i++)
        {
            memcpy(pabyData + 4 + i * 24, &(paoPoints[i].x), 8);
            memcpy(pabyData + 4 + i * 24 + 8, &(paoPoints[i].y), 8);
            if (padfZ == nullptr)
                memset(pabyData + 4 + i * 24 + 16, 0, 8);
            else
                memcpy(pabyData + 4 + i * 24 + 16, padfZ + i, 8);
        }
        OGRRoundCoordinatesIEEE754XYValues<24>(
            psOptions->sPrecision.nXYBitPrecision, pabyData + 4, nPointCount);
        OGRRoundCoordinatesIEEE754<24>(psOptions->sPrecision.nZBitPrecision,
                                       pabyData + 4 + 2 * sizeof(uint64_t),
                                       nPointCount);
    }
    else
    {
        nWords = 2 * static_cast<size_t>(nPointCount);
        memcpy(pabyData + 4, paoPoints, 16 * static_cast<size_t>(nPointCount));
        OGRRoundCoordinatesIEEE754XYValues<16>(
            psOptions->sPrecision.nXYBitPrecision, pabyData + 4, nPointCount);
    }

    /* -------------------------------------------------------------------- */
    /*      Swap if needed.                                                 */
    /* -------------------------------------------------------------------- */
    if (OGR_SWAP(psOptions->eByteOrder))
    {
        const int nCount = CPL_SWAP32(nPointCount);
        memcpy(pabyData, &nCount, 4);

        for (size_t i = 0; i < nWords; i++)
        {
            CPL_SWAPDOUBLE(pabyData + 4 + 8 * i);
        }
    }

    return OGRERR_NONE;
}

/************************************************************************/
/*                              _WkbSize()                              */
/*                                                                      */
/*      Helper method for OGRPolygon.  NOT THE NORMAL WkbSize() METHOD. */
/************************************************************************/

size_t OGRLinearRing::_WkbSize(int _flags) const

{
    if ((_flags & OGR_G_3D) && (_flags & OGR_G_MEASURED))
        return 4 + 32 * static_cast<size_t>(nPointCount);
    else if ((_flags & OGR_G_3D) || (_flags & OGR_G_MEASURED))
        return 4 + 24 * static_cast<size_t>(nPointCount);
    else
        return 4 + 16 * static_cast<size_t>(nPointCount);
}

//! @endcond

/************************************************************************/
/*                               clone()                                */
/*                                                                      */
/*      We override the OGRCurve clone() to ensure that we get the      */
/*      correct virtual table.                                          */
/************************************************************************/

OGRLinearRing *OGRLinearRing::clone() const

{
    OGRLinearRing *poNewLinearRing = new OGRLinearRing();
    poNewLinearRing->assignSpatialReference(getSpatialReference());

    poNewLinearRing->setPoints(nPointCount, paoPoints, padfZ, padfM);
    poNewLinearRing->flags = flags;

    return poNewLinearRing;
}

/************************************************************************/
/*                             reverseWindingOrder()                    */
/************************************************************************/

//! @cond Doxygen_Suppress
/** Reverse order of points.
 */
void OGRLinearRing::reverseWindingOrder()

{
    reversePoints();
}

//! @endcond

/************************************************************************/
/*                             closeRing()                              */
/************************************************************************/

void OGRLinearRing::closeRings()

{
    if (nPointCount < 2)
        return;

    if (getX(0) != getX(nPointCount - 1) || getY(0) != getY(nPointCount - 1) ||
        getZ(0) != getZ(nPointCount - 1))
    {
        OGRPoint oFirstPoint;
        getPoint(0, &oFirstPoint);
        addPoint(&oFirstPoint);
    }
}

/************************************************************************/
/*                              isPointInRing()                         */
/************************************************************************/

/** Returns whether the point is inside the ring.
 * @param poPoint point
 * @param bTestEnvelope set to TRUE if the presence of the point inside the
 *                      ring envelope must be checked first.
 * @return TRUE or FALSE.
 */
OGRBoolean OGRLinearRing::isPointInRing(const OGRPoint *poPoint,
                                        int bTestEnvelope) const
{
    if (nullptr == poPoint)
    {
        CPLDebug("OGR",
                 "OGRLinearRing::isPointInRing(const OGRPoint* poPoint) - "
                 "passed point is NULL!");
        return FALSE;
    }
    if (poPoint->IsEmpty())
    {
        return FALSE;
    }

    const int iNumPoints = getNumPoints();

    // Simple validation
    if (iNumPoints < 4)
        return FALSE;

    const double dfTestX = poPoint->getX();
    const double dfTestY = poPoint->getY();

    // Fast test if point is inside extent of the ring.
    if (bTestEnvelope)
    {
        OGREnvelope extent;
        getEnvelope(&extent);
        if (!(dfTestX >= extent.MinX && dfTestX <= extent.MaxX &&
              dfTestY >= extent.MinY && dfTestY <= extent.MaxY))
        {
            return FALSE;
        }
    }

    // For every point p in ring,
    // test if ray starting from given point crosses segment (p - 1, p)
    int iNumCrossings = 0;

    double prev_diff_x = getX(0) - dfTestX;
    double prev_diff_y = getY(0) - dfTestY;

    for (int iPoint = 1; iPoint < iNumPoints; iPoint++)
    {
        const double x1 = getX(iPoint) - dfTestX;
        const double y1 = getY(iPoint) - dfTestY;

        const double x2 = prev_diff_x;
        const double y2 = prev_diff_y;

        if (((y1 > 0) && (y2 <= 0)) || ((y2 > 0) && (y1 <= 0)))
        {
            // Check if ray intersects with segment of the ring
            const double dfIntersection = (x1 * y2 - x2 * y1) / (y2 - y1);
            if (0.0 < dfIntersection)
            {
                // Count intersections
                iNumCrossings++;
            }
        }

        prev_diff_x = x1;
        prev_diff_y = y1;
    }

    // If iNumCrossings number is even, given point is outside the ring,
    // when the crossings number is odd, the point is inside the ring.
    return iNumCrossings % 2;  // OGRBoolean
}

/************************************************************************/
/*                       isPointOnRingBoundary()                        */
/************************************************************************/

/** Returns whether the point is on the ring boundary.
 * @param poPoint point
 * @param bTestEnvelope set to TRUE if the presence of the point inside the
 *                      ring envelope must be checked first.
 * @return TRUE or FALSE.
 */
OGRBoolean OGRLinearRing::isPointOnRingBoundary(const OGRPoint *poPoint,
                                                int bTestEnvelope) const
{
    if (nullptr == poPoint)
    {
        CPLDebug("OGR", "OGRLinearRing::isPointOnRingBoundary(const OGRPoint* "
                        "poPoint) - passed point is NULL!");
        return 0;
    }

    const int iNumPoints = getNumPoints();

    // Simple validation.
    if (iNumPoints < 4)
        return 0;

    const double dfTestX = poPoint->getX();
    const double dfTestY = poPoint->getY();

    // Fast test if point is inside extent of the ring
    if (bTestEnvelope)
    {
        OGREnvelope extent;
        getEnvelope(&extent);
        if (!(dfTestX >= extent.MinX && dfTestX <= extent.MaxX &&
              dfTestY >= extent.MinY && dfTestY <= extent.MaxY))
        {
            return 0;
        }
    }

    double prev_diff_x = dfTestX - getX(0);
    double prev_diff_y = dfTestY - getY(0);

    for (int iPoint = 1; iPoint < iNumPoints; iPoint++)
    {
        const double dx1 = dfTestX - getX(iPoint);
        const double dy1 = dfTestY - getY(iPoint);

        const double dx2 = prev_diff_x;
        const double dy2 = prev_diff_y;

        // If the point is on the segment, return immediately.
        // FIXME? If the test point is not exactly identical to one of
        // the vertices of the ring, but somewhere on a segment, there's
        // little chance that we get 0. So that should be tested against some
        // epsilon.

        if (dx1 * dy2 - dx2 * dy1 == 0)
        {
            // If iPoint and iPointPrev are the same, go on.
            if (!(dx1 == dx2 && dy1 == dy2))
            {
                const double dx_segment = getX(iPoint) - getX(iPoint - 1);
                const double dy_segment = getY(iPoint) - getY(iPoint - 1);
                const double crossproduct = dx2 * dx_segment + dy2 * dy_segment;
                if (crossproduct >= 0)
                {
                    const double sq_length_seg =
                        dx_segment * dx_segment + dy_segment * dy_segment;
                    if (crossproduct <= sq_length_seg)
                    {
                        return 1;
                    }
                }
            }
        }

        prev_diff_x = dx1;
        prev_diff_y = dy1;
    }

    return 0;
}

/************************************************************************/
/*                             transform()                              */
/************************************************************************/

OGRErr OGRLinearRing::transform(OGRCoordinateTransformation *poCT)

{
    const bool bIsClosed = getNumPoints() > 2 && CPL_TO_BOOL(get_IsClosed());
    OGRErr eErr = OGRLineString::transform(poCT);
    if (bIsClosed && eErr == OGRERR_NONE && !get_IsClosed())
    {
        CPLDebug("OGR", "Linearring is not closed after coordinate "
                        "transformation. Forcing last point to be identical to "
                        "first one");
        // Force last point to be identical to first point.
        // This is a safety belt in case the reprojection of the same coordinate
        // isn't perfectly stable. This can for example happen in very rare
        // cases when reprojecting a cutline with a RPC transform with a DEM
        // that is a VRT whose sources are resampled...
        OGRPoint oStartPoint;
        StartPoint(&oStartPoint);

        setPoint(getNumPoints() - 1, &oStartPoint);
    }
    return eErr;
}

/************************************************************************/
/*                          CastToLineString()                          */
/************************************************************************/

/**
 * \brief Cast to line string.
 *
 * The passed in geometry is consumed and a new one returned .
 *
 * @param poLR the input geometry - ownership is passed to the method.
 * @return new geometry.
 */

OGRLineString *OGRLinearRing::CastToLineString(OGRLinearRing *poLR)
{
    return TransferMembersAndDestroy(poLR, new OGRLineString());
}

//! @cond Doxygen_Suppress
/************************************************************************/
/*                     GetCasterToLineString()                          */
/************************************************************************/

OGRLineString *OGRLinearRing::CasterToLineString(OGRCurve *poCurve)
{
    return OGRLinearRing::CastToLineString(poCurve->toLinearRing());
}

OGRCurveCasterToLineString OGRLinearRing::GetCasterToLineString() const
{
    return OGRLinearRing::CasterToLineString;
}

/************************************************************************/
/*                        GetCasterToLinearRing()                       */
/************************************************************************/

static OGRLinearRing *CasterToLinearRing(OGRCurve *poCurve)
{
    return poCurve->toLinearRing();
}

OGRCurveCasterToLinearRing OGRLinearRing::GetCasterToLinearRing() const
{
    return ::CasterToLinearRing;
}

//! @endcond

Coverage Report

Created: 2025-06-13 06:18