qgslinestring.cpp

Go to the documentation of this file.

/***************************************************************************
                         qgslinestring.cpp
                         -------------------
    begin                : September 2014
    copyright            : (C) 2014 by Marco Hugentobler
    email                : marco at sourcepole dot ch
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
 
#include "qgslinestring.h"
 
#include <cmath>
#include <limits>
#include <memory>
#include <nlohmann/json.hpp>
 
#include "qgsbox3d.h"
#include "qgscompoundcurve.h"
#include "qgscoordinatetransform.h"
#include "qgsgeometryutils.h"
#include "qgsgeometryutils_base.h"
#include "qgslinesegment.h"
#include "qgsvector3d.h"
 
#include <QDomDocument>
#include <QJsonObject>
#include <QPainter>
#include <QString>
 
using namespace Qt::StringLiterals;
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
QgsLineString::QgsLineString()
{
  mWkbType = Qgis::WkbType::LineString;
}
 
QgsLineString::QgsLineString( const QVector<QgsPoint> &points )
{
  if ( points.isEmpty() )
  {
    mWkbType = Qgis::WkbType::LineString;
    return;
  }
  Qgis::WkbType ptType = points.at( 0 ).wkbType();
  mWkbType = QgsWkbTypes::zmType( Qgis::WkbType::LineString, QgsWkbTypes::hasZ( ptType ), QgsWkbTypes::hasM( ptType ) );
  mX.resize( points.count() );
  mY.resize( points.count() );
  double *x = mX.data();
  double *y = mY.data();
  double *z = nullptr;
  double *m = nullptr;
  if ( QgsWkbTypes::hasZ( mWkbType ) )
  {
    mZ.resize( points.count() );
    z = mZ.data();
  }
  if ( QgsWkbTypes::hasM( mWkbType ) )
  {
    mM.resize( points.count() );
    m = mM.data();
  }
 
  for ( const QgsPoint &pt : points )
  {
    *x++ = pt.x();
    *y++ = pt.y();
    if ( z )
      *z++ = pt.z();
    if ( m )
      *m++ = pt.m();
  }
}
 
QgsLineString::QgsLineString( const QVector<double> &x, const QVector<double> &y, const QVector<double> &z, const QVector<double> &m, bool is25DType )
{
  mWkbType = Qgis::WkbType::LineString;
  int pointCount = std::min( x.size(), y.size() );
  if ( x.size() == pointCount )
  {
    mX = x;
  }
  else
  {
    mX = x.mid( 0, pointCount );
  }
  if ( y.size() == pointCount )
  {
    mY = y;
  }
  else
  {
    mY = y.mid( 0, pointCount );
  }
  if ( !z.isEmpty() && z.count() >= pointCount )
  {
    mWkbType = is25DType ? Qgis::WkbType::LineString25D : Qgis::WkbType::LineStringZ;
    if ( z.size() == pointCount )
    {
      mZ = z;
    }
    else
    {
      mZ = z.mid( 0, pointCount );
    }
  }
  if ( !m.isEmpty() && m.count() >= pointCount )
  {
    mWkbType = QgsWkbTypes::addM( mWkbType );
    if ( m.size() == pointCount )
    {
      mM = m;
    }
    else
    {
      mM = m.mid( 0, pointCount );
    }
  }
}
 
QgsLineString::QgsLineString( const QgsPoint &p1, const QgsPoint &p2 )
{
  mWkbType = Qgis::WkbType::LineString;
  mX.resize( 2 );
  mX[0] = p1.x();
  mX[1] = p2.x();
  mY.resize( 2 );
  mY[0] = p1.y();
  mY[1] = p2.y();
  if ( p1.is3D() )
  {
    mWkbType = QgsWkbTypes::addZ( mWkbType );
    mZ.resize( 2 );
    mZ[0] = p1.z();
    mZ[1] = p2.z();
  }
  if ( p1.isMeasure() )
  {
    mWkbType = QgsWkbTypes::addM( mWkbType );
    mM.resize( 2 );
    mM[0] = p1.m();
    mM[1] = p2.m();
  }
}
 
QgsLineString::QgsLineString( const QVector<QgsPointXY> &points )
{
  mWkbType = Qgis::WkbType::LineString;
  mX.reserve( points.size() );
  mY.reserve( points.size() );
  for ( const QgsPointXY &p : points )
  {
    mX << p.x();
    mY << p.y();
  }
}
 
QgsLineString::QgsLineString( const QgsLineSegment2D &segment )
{
  mWkbType = Qgis::WkbType::LineString;
  mX.resize( 2 );
  mY.resize( 2 );
  mX[0] = segment.startX();
  mX[1] = segment.endX();
  mY[0] = segment.startY();
  mY[1] = segment.endY();
}
 
std::unique_ptr< QgsLineString > QgsLineString::fromBezierCurve( const QgsPoint &start, const QgsPoint &controlPoint1, const QgsPoint &controlPoint2, const QgsPoint &end, int segments )
{
  if ( segments == 0 )
    return std::make_unique< QgsLineString >();
 
  QVector<double> x( segments + 1 );
  QVector<double> y( segments + 1 );
  QVector<double> z;
  QVector<double> m;
 
  const bool hasZ = start.is3D() && controlPoint1.is3D() && controlPoint2.is3D() && end.is3D();
  if ( hasZ )
  {
    z.resize( segments + 1 );
  }
 
  const bool hasM = start.isMeasure() && controlPoint1.isMeasure() && controlPoint2.isMeasure() && end.isMeasure();
  if ( hasM )
  {
    m.resize( segments + 1 );
  }
 
  double *xData = x.data();
  double *yData = y.data();
  double *zData = z.data(); // will be nullptr if !hasZ
  double *mData = m.data(); // will be nullptr if !hasM
 
  const double step = 1.0 / segments;
 
  for ( int i = 0; i <= segments; ++i )
  {
    const double t = i * step;
 
    double ix, iy; // interpolated x, y
    double iz = std::numeric_limits<double>::quiet_NaN();
    double im = std::numeric_limits<double>::quiet_NaN();
 
    QgsGeometryUtilsBase::interpolatePointOnCubicBezier(
      start.x(),
      start.y(),
      start.z(),
      start.m(),
      controlPoint1.x(),
      controlPoint1.y(),
      controlPoint1.z(),
      controlPoint1.m(),
      controlPoint2.x(),
      controlPoint2.y(),
      controlPoint2.z(),
      controlPoint2.m(),
      end.x(),
      end.y(),
      end.z(),
      end.m(),
      t,
      hasZ,
      hasM,
      ix,
      iy,
      iz,
      im
    );
 
    *xData++ = ix;
    *yData++ = iy;
    if ( hasZ )
      *zData++ = iz;
    if ( hasM )
      *mData++ = im;
  }
 
  return std::make_unique< QgsLineString >( x, y, z, m );
}
 
std::unique_ptr< QgsLineString > QgsLineString::fromQPolygonF( const QPolygonF &polygon )
{
  QVector< double > x;
  QVector< double > y;
  x.resize( polygon.count() );
  y.resize( polygon.count() );
  double *xData = x.data();
  double *yData = y.data();
 
  const QPointF *src = polygon.data();
  for ( int i = 0; i < polygon.size(); ++i )
  {
    *xData++ = src->x();
    *yData++ = src->y();
    src++;
  }
 
  return std::make_unique< QgsLineString >( x, y );
}
 
QgsLineString *QgsLineString::clone() const
{
  return new QgsLineString( *this );
}
 
void QgsLineString::clear()
{
  mWkbType = Qgis::WkbType::LineString;
  QgsSimpleCurve::clear();
}
 
int QgsLineString::indexOf( const QgsPoint &point ) const
{
  const int size = mX.size();
  if ( size == 0 )
    return -1;
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  const bool useZ = is3D();
  const bool useM = isMeasure();
  const double *z = useZ ? mZ.constData() : nullptr;
  const double *m = useM ? mM.constData() : nullptr;
 
  for ( int i = 0; i < size; ++i )
  {
    if ( qgsDoubleNear( *x, point.x() ) && qgsDoubleNear( *y, point.y() ) && ( !useZ || qgsDoubleNear( *z, point.z() ) ) && ( !useM || qgsDoubleNear( *m, point.m() ) ) )
      return i;
 
    x++;
    y++;
    if ( useZ )
      z++;
    if ( useM )
      m++;
  }
  return -1;
}
 
bool QgsLineString::isValid( QString &error, Qgis::GeometryValidityFlags flags ) const
{
  if ( !isEmpty() && ( numPoints() < 2 ) )
  {
    error = QObject::tr( "LineString has less than 2 points and is not empty." );
    return false;
  }
  return QgsCurve::isValid( error, flags );
}
 
QgsLineString *QgsLineString::snappedToGrid( double hSpacing, double vSpacing, double dSpacing, double mSpacing, bool removeRedundantPoints ) const
{
  // prepare result
  std::unique_ptr<QgsLineString> result { createEmptyWithSameType() };
 
  bool res = snapToGridPrivate( hSpacing, vSpacing, dSpacing, mSpacing, mX, mY, mZ, mM, result->mX, result->mY, result->mZ, result->mM, removeRedundantPoints );
  if ( res )
    return result.release();
  else
    return nullptr;
}
 
bool QgsLineString::removeDuplicateNodes( double epsilon, bool useZValues )
{
  if ( mX.count() <= 2 )
    return false; // don't create degenerate lines
  bool result = false;
  double prevX = mX.at( 0 );
  double prevY = mY.at( 0 );
  bool hasZ = is3D();
  bool useZ = hasZ && useZValues;
  double prevZ = useZ ? mZ.at( 0 ) : 0;
  int i = 1;
  int remaining = mX.count();
  while ( i < remaining )
  {
    double currentX = mX.at( i );
    double currentY = mY.at( i );
    double currentZ = useZ ? mZ.at( i ) : 0;
    if ( qgsDoubleNear( currentX, prevX, epsilon ) && qgsDoubleNear( currentY, prevY, epsilon ) && ( !useZ || qgsDoubleNear( currentZ, prevZ, epsilon ) ) )
    {
      result = true;
      // remove point
      mX.removeAt( i );
      mY.removeAt( i );
      if ( hasZ )
        mZ.removeAt( i );
      remaining--;
    }
    else
    {
      prevX = currentX;
      prevY = currentY;
      prevZ = currentZ;
      i++;
    }
  }
  return result;
}
 
bool QgsLineString::isClosed2D() const
{
  if ( mX.empty() )
    return false;
 
  return qgsDoubleNear( mX.first(), mX.last() ) && qgsDoubleNear( mY.first(), mY.last() );
}
 
bool QgsLineString::isClosed() const
{
  bool closed = isClosed2D();
 
  if ( is3D() && closed )
    closed &= qgsDoubleNear( mZ.first(), mZ.last() ) || ( std::isnan( mZ.first() ) && std::isnan( mZ.last() ) );
  return closed;
}
 
// As `bool boundingBoxIntersects( const QgsBox3D &box3d )` and `bool boundingBoxIntersects( const QgsRectangle &rectangle )` are nearly
// the same: if one of these functions is changed then remember to also update the other accordingly
bool QgsLineString::boundingBoxIntersects( const QgsRectangle &rectangle ) const
{
  if ( mX.empty() )
    return false;
 
  if ( !mBoundingBox.isNull() )
  {
    return mBoundingBox.toRectangle().intersects( rectangle );
  }
  const int nb = mX.size();
 
  // We are a little fancy here!
  if ( nb > 40 )
  {
    // if a large number of vertices, take some sample vertices at 1/5th increments through the linestring
    // and test whether any are inside the rectangle. Maybe we can shortcut a lot of iterations by doing this!
    // (why 1/5th? it's picked so that it works nicely for polygon rings which are almost rectangles, so the vertex extremities
    // will fall on approximately these vertex indices)
    if ( rectangle.contains( mX.at( 0 ), mY.at( 0 ) )
         || rectangle.contains( mX.at( static_cast< int >( nb * 0.2 ) ), mY.at( static_cast< int >( nb * 0.2 ) ) )
         || rectangle.contains( mX.at( static_cast< int >( nb * 0.4 ) ), mY.at( static_cast< int >( nb * 0.4 ) ) )
         || rectangle.contains( mX.at( static_cast< int >( nb * 0.6 ) ), mY.at( static_cast< int >( nb * 0.6 ) ) )
         || rectangle.contains( mX.at( static_cast< int >( nb * 0.8 ) ), mY.at( static_cast< int >( nb * 0.8 ) ) )
         || rectangle.contains( mX.at( nb - 1 ), mY.at( nb - 1 ) ) )
      return true;
  }
 
  // Be even MORE fancy! Given that bounding box calculation is non-free, cached, and we don't
  // already have it, we start performing the bounding box calculation while we are testing whether
  // each point falls inside the rectangle. That way if we end up testing the majority of the points
  // anyway, we can update the cached bounding box with the results we've calculated along the way
  // and save future calls to calculate the bounding box!
  double xmin = std::numeric_limits<double>::max();
  double ymin = std::numeric_limits<double>::max();
  double zmin = -std::numeric_limits<double>::max();
  double xmax = -std::numeric_limits<double>::max();
  double ymax = -std::numeric_limits<double>::max();
  double zmax = -std::numeric_limits<double>::max();
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  const double *z = is3D() ? mZ.constData() : nullptr;
  bool foundPointInRectangle = false;
  for ( int i = 0; i < nb; ++i )
  {
    const double px = *x++;
    xmin = std::min( xmin, px );
    xmax = std::max( xmax, px );
    const double py = *y++;
    ymin = std::min( ymin, py );
    ymax = std::max( ymax, py );
    if ( z )
    {
      const double pz = *z++;
      zmin = std::min( zmin, pz );
      zmax = std::max( zmax, pz );
    }
 
    if ( !foundPointInRectangle && rectangle.contains( px, py ) )
    {
      foundPointInRectangle = true;
 
      // now... we have a choice to make. If we've already looped through the majority of the points
      // in this linestring then let's just continue to iterate through the remainder so that we can
      // complete the overall bounding box calculation we've already mostly done. If however we're only
      // just at the start of iterating the vertices, we shortcut out early and leave the bounding box
      // uncalculated
      if ( i < nb * 0.5 )
        return true;
    }
  }
 
  // at this stage we now know the overall bounding box of the linestring, so let's cache
  // it so we don't ever have to calculate this again. We've done all the hard work anyway!
  mBoundingBox = QgsBox3D( xmin, ymin, zmin, xmax, ymax, zmax, false );
 
  if ( foundPointInRectangle )
    return true;
 
  // NOTE: if none of the points in the line actually fell inside the rectangle, it doesn't
  // exclude that the OVERALL bounding box of the linestring itself intersects the rectangle!!
  // So we fall back to the parent class method which compares the overall bounding box against
  // the rectangle... and this will be very cheap now that we've already calculated and cached
  // the linestring's bounding box!
  return QgsCurve::boundingBoxIntersects( rectangle );
}
 
// As `bool boundingBoxIntersects( const QgsBox3D &box3d )` and `bool boundingBoxIntersects( const QgsRectangle &rectangle )` are nearly
// the same: if one of these functions is changed then remember to also update the other accordingly
bool QgsLineString::boundingBoxIntersects( const QgsBox3D &box3d ) const
{
  if ( mX.empty() )
    return false;
 
  if ( mZ.empty() )
    return boundingBoxIntersects( box3d.toRectangle() );
 
  if ( !mBoundingBox.isNull() )
  {
    return mBoundingBox.intersects( box3d );
  }
  const int nb = mX.size();
 
  // We are a little fancy here!
  if ( nb > 40 )
  {
    // if a large number of vertices, take some sample vertices at 1/5th increments through the linestring
    // and test whether any are inside the rectangle. Maybe we can shortcut a lot of iterations by doing this!
    // (why 1/5th? it's picked so that it works nicely for polygon rings which are almost rectangles, so the vertex extremities
    // will fall on approximately these vertex indices)
    if ( box3d.contains( mX.at( 0 ), mY.at( 0 ), mZ.at( 0 ) )
         || box3d.contains( mX.at( static_cast< int >( nb * 0.2 ) ), mY.at( static_cast< int >( nb * 0.2 ) ), mZ.at( static_cast< int >( nb * 0.2 ) ) )
         || box3d.contains( mX.at( static_cast< int >( nb * 0.4 ) ), mY.at( static_cast< int >( nb * 0.4 ) ), mZ.at( static_cast< int >( nb * 0.4 ) ) )
         || box3d.contains( mX.at( static_cast< int >( nb * 0.6 ) ), mY.at( static_cast< int >( nb * 0.6 ) ), mZ.at( static_cast< int >( nb * 0.6 ) ) )
         || box3d.contains( mX.at( static_cast< int >( nb * 0.8 ) ), mY.at( static_cast< int >( nb * 0.8 ) ), mZ.at( static_cast< int >( nb * 0.8 ) ) )
         || box3d.contains( mX.at( nb - 1 ), mY.at( nb - 1 ), mZ.at( nb - 1 ) ) )
      return true;
  }
 
  // Be even MORE fancy! Given that bounding box calculation is non-free, cached, and we don't
  // already have it, we start performing the bounding box calculation while we are testing whether
  // each point falls inside the rectangle. That way if we end up testing the majority of the points
  // anyway, we can update the cached bounding box with the results we've calculated along the way
  // and save future calls to calculate the bounding box!
  double xmin = std::numeric_limits<double>::max();
  double ymin = std::numeric_limits<double>::max();
  double zmin = std::numeric_limits<double>::max();
  double xmax = -std::numeric_limits<double>::max();
  double ymax = -std::numeric_limits<double>::max();
  double zmax = -std::numeric_limits<double>::max();
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  const double *z = mZ.constData();
  bool foundPointInBox = false;
  for ( int i = 0; i < nb; ++i )
  {
    const double px = *x++;
    xmin = std::min( xmin, px );
    xmax = std::max( xmax, px );
    const double py = *y++;
    ymin = std::min( ymin, py );
    ymax = std::max( ymax, py );
    const double pz = *z++;
    zmin = std::min( zmin, pz );
    zmax = std::max( zmax, pz );
 
    if ( !foundPointInBox && box3d.contains( px, py, pz ) )
    {
      foundPointInBox = true;
 
      // now... we have a choice to make. If we've already looped through the majority of the points
      // in this linestring then let's just continue to iterate through the remainder so that we can
      // complete the overall bounding box calculation we've already mostly done. If however we're only
      // just at the start of iterating the vertices, we shortcut out early and leave the bounding box
      // uncalculated
      if ( i < nb * 0.5 )
        return true;
    }
  }
 
  // at this stage we now know the overall bounding box of the linestring, so let's cache
  // it so we don't ever have to calculate this again. We've done all the hard work anyway!
  mBoundingBox = QgsBox3D( xmin, ymin, zmin, xmax, ymax, zmax, false );
 
  if ( foundPointInBox )
    return true;
 
  // NOTE: if none of the points in the line actually fell inside the rectangle, it doesn't
  // exclude that the OVERALL bounding box of the linestring itself intersects the rectangle!!
  // So we fall back to the parent class method which compares the overall bounding box against
  // the rectangle... and this will be very cheap now that we've already calculated and cached
  // the linestring's bounding box!
  return QgsCurve::boundingBoxIntersects( box3d );
}
 
QVector< QgsVertexId > QgsLineString::collectDuplicateNodes( double epsilon, bool useZValues ) const
{
  QVector< QgsVertexId > res;
  if ( mX.count() <= 1 )
    return res;
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  bool hasZ = is3D();
  bool useZ = hasZ && useZValues;
  const double *z = useZ ? mZ.constData() : nullptr;
 
  double prevX = *x++;
  double prevY = *y++;
  double prevZ = z ? *z++ : 0;
 
  QgsVertexId id;
  for ( int i = 1; i < mX.count(); ++i )
  {
    double currentX = *x++;
    double currentY = *y++;
    double currentZ = useZ ? *z++ : 0;
    if ( qgsDoubleNear( currentX, prevX, epsilon ) && qgsDoubleNear( currentY, prevY, epsilon ) && ( !useZ || qgsDoubleNear( currentZ, prevZ, epsilon ) ) )
    {
      id.vertex = i;
      res << id;
    }
    else
    {
      prevX = currentX;
      prevY = currentY;
      prevZ = currentZ;
    }
  }
  return res;
}
 
QPolygonF QgsLineString::asQPolygonF() const
{
  const int nb = mX.size();
  QPolygonF points( nb );
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  QPointF *dest = points.data();
  for ( int i = 0; i < nb; ++i )
  {
    *dest++ = QPointF( *x++, *y++ );
  }
  return points;
}
 
 
void simplifySection( int i, int j, const double *x, const double *y, std::vector< bool > &usePoint, const double distanceToleranceSquared, const double epsilon )
{
  if ( i + 1 == j )
  {
    return;
  }
 
  double maxDistanceSquared = -1.0;
 
  int maxIndex = i;
  double mx, my;
 
  for ( int k = i + 1; k < j; k++ )
  {
    const double distanceSquared = QgsGeometryUtilsBase::sqrDistToLine( x[k], y[k], x[i], y[i], x[j], y[j], mx, my, epsilon );
 
    if ( distanceSquared > maxDistanceSquared )
    {
      maxDistanceSquared = distanceSquared;
      maxIndex = k;
    }
  }
  if ( maxDistanceSquared <= distanceToleranceSquared )
  {
    for ( int k = i + 1; k < j; k++ )
    {
      usePoint[k] = false;
    }
  }
  else
  {
    simplifySection( i, maxIndex, x, y, usePoint, distanceToleranceSquared, epsilon );
    simplifySection( maxIndex, j, x, y, usePoint, distanceToleranceSquared, epsilon );
  }
};
 
QgsLineString *QgsLineString::simplifyByDistance( double tolerance ) const
{
  if ( mX.empty() )
  {
    return new QgsLineString();
  }
 
  // ported from GEOS DouglasPeuckerLineSimplifier::simplify
 
  const double distanceToleranceSquared = tolerance * tolerance;
  const double *xData = mX.constData();
  const double *yData = mY.constData();
  const double *zData = mZ.constData();
  const double *mData = mM.constData();
 
  const int size = mX.size();
 
  std::vector< bool > usePoint( size, true );
 
  constexpr double epsilon = 4 * std::numeric_limits<double>::epsilon();
  simplifySection( 0, size - 1, xData, yData, usePoint, distanceToleranceSquared, epsilon );
 
  QVector< double > newX;
  newX.reserve( size );
  QVector< double > newY;
  newY.reserve( size );
 
  const bool hasZ = is3D();
  const bool hasM = isMeasure();
  QVector< double > newZ;
  if ( hasZ )
    newZ.reserve( size );
  QVector< double > newM;
  if ( hasM )
    newM.reserve( size );
 
  for ( int i = 0, n = size; i < n; ++i )
  {
    if ( usePoint[i] || i == n - 1 )
    {
      newX.append( xData[i] );
      newY.append( yData[i] );
      if ( hasZ )
        newZ.append( zData[i] );
      if ( hasM )
        newM.append( mData[i] );
    }
  }
 
  const bool simplifyRing = isRing();
  const int newSize = newX.size();
  if ( simplifyRing && newSize > 3 )
  {
    double mx, my;
    const double distanceSquared = QgsGeometryUtilsBase::sqrDistToLine( newX[0], newY[0], newX[newSize - 2], newY[newSize - 2], newX[1], newY[1], mx, my, epsilon );
 
    if ( distanceSquared <= distanceToleranceSquared )
    {
      newX.removeFirst();
      newX.last() = newX.first();
      newY.removeFirst();
      newY.last() = newY.first();
      if ( hasZ )
      {
        newZ.removeFirst();
        newZ.last() = newZ.first();
      }
      if ( hasM )
      {
        newM.removeFirst();
        newM.last() = newM.first();
      }
    }
  }
 
  return new QgsLineString( newX, newY, newZ, newM );
}
 
QgsBox3D QgsLineString::calculateBoundingBox3D() const
{
  if ( mX.empty() )
  {
    return QgsBox3D();
  }
 
  auto result2D = std::minmax_element( mX.begin(), mX.end() );
  const double xmin = *result2D.first;
  const double xmax = *result2D.second;
  result2D = std::minmax_element( mY.begin(), mY.end() );
  const double ymin = *result2D.first;
  const double ymax = *result2D.second;
 
  double zmin = std::numeric_limits< double >::quiet_NaN();
  double zmax = std::numeric_limits< double >::quiet_NaN();
 
  if ( is3D() )
  {
    auto resultZ = std::minmax_element( mZ.begin(), mZ.end() );
    zmin = *resultZ.first;
    zmax = *resultZ.second;
  }
 
  return QgsBox3D( xmin, ymin, zmin, xmax, ymax, zmax );
}
 
QgsBox3D QgsLineString::calculateBoundingBox3d() const
{
  return calculateBoundingBox3D();
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
QDomElement QgsLineString::asGml2( QDomDocument &doc, int precision, const QString &ns, const AxisOrder axisOrder ) const
{
  QgsPointSequence pts;
  points( pts );
 
  QDomElement elemLineString = doc.createElementNS( ns, u"LineString"_s );
 
  if ( isEmpty() )
    return elemLineString;
 
  elemLineString.appendChild( QgsGeometryUtils::pointsToGML2( pts, doc, precision, ns, axisOrder ) );
 
  return elemLineString;
}
 
QDomElement QgsLineString::asGml3( QDomDocument &doc, int precision, const QString &ns, const QgsAbstractGeometry::AxisOrder axisOrder ) const
{
  QgsPointSequence pts;
  points( pts );
 
  QDomElement elemLineString = doc.createElementNS( ns, u"LineString"_s );
 
  if ( isEmpty() )
    return elemLineString;
 
  elemLineString.appendChild( QgsGeometryUtils::pointsToGML3( pts, doc, precision, ns, is3D(), axisOrder ) );
  return elemLineString;
}
 
json QgsLineString::asJsonObject( int precision ) const
{
  QgsPointSequence pts;
  points( pts );
  return { { "type", "LineString" }, { "coordinates", QgsGeometryUtils::pointsToJson( pts, precision ) } };
}
 
QString QgsLineString::asKml( int precision ) const
{
  QString kml;
  if ( isRing() )
  {
    kml.append( "<LinearRing>"_L1 );
  }
  else
  {
    kml.append( "<LineString>"_L1 );
  }
  bool z = is3D();
  kml.append( "<altitudeMode>"_L1 );
  if ( z )
  {
    kml.append( "absolute"_L1 );
  }
  else
  {
    kml.append( "clampToGround"_L1 );
  }
  kml.append( "</altitudeMode>"_L1 );
  kml.append( "<coordinates>"_L1 );
 
  int nPoints = mX.size();
  for ( int i = 0; i < nPoints; ++i )
  {
    if ( i > 0 )
    {
      kml.append( " "_L1 );
    }
    kml.append( qgsDoubleToString( mX[i], precision ) );
    kml.append( ","_L1 );
    kml.append( qgsDoubleToString( mY[i], precision ) );
    if ( z )
    {
      kml.append( ","_L1 );
      kml.append( qgsDoubleToString( mZ[i], precision ) );
    }
    else
    {
      kml.append( ",0"_L1 );
    }
  }
  kml.append( "</coordinates>"_L1 );
  if ( isRing() )
  {
    kml.append( "</LinearRing>"_L1 );
  }
  else
  {
    kml.append( "</LineString>"_L1 );
  }
  return kml;
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
double QgsLineString::length() const
{
  double total = 0;
  const int size = mX.size();
  if ( size < 2 )
    return 0;
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  double dx, dy;
 
  double prevX = *x++;
  double prevY = *y++;
 
  for ( int i = 1; i < size; ++i )
  {
    dx = *x - prevX;
    dy = *y - prevY;
    total += std::sqrt( dx * dx + dy * dy );
 
    prevX = *x++;
    prevY = *y++;
  }
  return total;
}
 
std::tuple<std::unique_ptr<QgsCurve>, std::unique_ptr<QgsCurve> > QgsLineString::splitCurveAtVertex( int index ) const
{
  QVector< double > x1, y1, z1, m1;
  QVector< double > x2, y2, z2, m2;
  QgsSimpleCurve::splitCurveAtVertexProtected( index, x1, y1, z1, m1, x2, y2, z2, m2 );
 
  std::unique_ptr< QgsLineString > first;
  if ( x1.isEmpty() || ( x1.size() < 2 && x2.size() >= 2 ) )
    first = std::make_unique< QgsLineString >();
  else
    first = std::make_unique< QgsLineString >( x1, y1, z1, m1 );
 
  std::unique_ptr< QgsLineString > second;
  if ( x2.isEmpty() || x2.size() < 2 )
    second = std::make_unique< QgsLineString >();
  else
    second = std::make_unique< QgsLineString >( x2, y2, z2, m2 );
 
  return std::make_tuple( std::move( first ), std::move( second ) );
}
 
QVector<QgsLineString *> QgsLineString::splitToDisjointXYParts() const
{
  const double *allPointsX = xData();
  const double *allPointsY = yData();
  size_t allPointsCount = numPoints();
  QVector<double> partX;
  QVector<double> partY;
  QSet<QgsPointXY> partPointSet;
 
  QVector<QgsLineString *> disjointParts;
  for ( size_t i = 0; i < allPointsCount; i++ )
  {
    const QgsPointXY point( *allPointsX++, *allPointsY++ );
    if ( partPointSet.contains( point ) )
    {
      // This point is used multiple times, cut the curve and add the
      // current part
      disjointParts.push_back( new QgsLineString( partX, partY ) );
      // Now start a new part containing the last line
      partX = { partX.last() };
      partY = { partY.last() };
      partPointSet = { QgsPointXY( partX[0], partY[0] ) };
    }
    partX.push_back( point.x() );
    partY.push_back( point.y() );
    partPointSet.insert( point );
  }
  // Add the last part (if we didn't stop by closing the loop)
  if ( partX.size() > 1 || disjointParts.size() == 0 )
    disjointParts.push_back( new QgsLineString( partX, partY ) );
 
  return disjointParts;
}
 
double QgsLineString::length3D() const
{
  if ( is3D() )
  {
    double total = 0;
    const int size = mX.size();
    if ( size < 2 )
      return 0;
 
    const double *x = mX.constData();
    const double *y = mY.constData();
    const double *z = mZ.constData();
    double dx, dy, dz;
 
    double prevX = *x++;
    double prevY = *y++;
    double prevZ = *z++;
 
    for ( int i = 1; i < size; ++i )
    {
      dx = *x - prevX;
      dy = *y - prevY;
      dz = *z - prevZ;
      total += std::sqrt( dx * dx + dy * dy + dz * dz );
 
      prevX = *x++;
      prevY = *y++;
      prevZ = *z++;
    }
    return total;
  }
  else
  {
    return length();
  }
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
QgsLineString *QgsLineString::curveToLine( double tolerance, SegmentationToleranceType toleranceType ) const
{
  Q_UNUSED( tolerance )
  Q_UNUSED( toleranceType )
  return clone();
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
void QgsLineString::setPoints( size_t size, const double *x, const double *y, const double *z, const double *m )
{
  clearCache(); //set bounding box invalid
 
  if ( size == 0 )
  {
    clear();
    return;
  }
 
  const bool hasZ = static_cast< bool >( z );
  const bool hasM = static_cast< bool >( m );
 
  if ( hasZ && hasM )
  {
    mWkbType = Qgis::WkbType::LineStringZM;
  }
  else if ( hasZ )
  {
    mWkbType = Qgis::WkbType::LineStringZ;
  }
  else if ( hasM )
  {
    mWkbType = Qgis::WkbType::LineStringM;
  }
  else
  {
    mWkbType = Qgis::WkbType::LineString;
  }
 
  mX.resize( size );
  mY.resize( size );
  double *destX = mX.data();
  double *destY = mY.data();
  double *destZ = nullptr;
  if ( hasZ )
  {
    mZ.resize( size );
    destZ = mZ.data();
  }
  else
  {
    mZ.clear();
  }
  double *destM = nullptr;
  if ( hasM )
  {
    mM.resize( size );
    destM = mM.data();
  }
  else
  {
    mM.clear();
  }
 
  for ( size_t i = 0; i < size; ++i )
  {
    *destX++ = *x++;
    *destY++ = *y++;
    if ( hasZ )
    {
      *destZ++ = *z++;
    }
    if ( hasM )
    {
      *destM++ = *m++;
    }
  }
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
QgsLineString *QgsLineString::reversed() const
{
  return qgis::down_cast< QgsLineString *>( QgsSimpleCurve::reversed() );
}
 
void QgsLineString::visitPointsByRegularDistance(
  const double distance, const std::function<bool( double, double, double, double, double, double, double, double, double, double, double, double )> &visitPoint
) const
{
  if ( distance < 0 )
    return;
 
  double distanceTraversed = 0;
  const int totalPoints = numPoints();
  if ( totalPoints == 0 )
    return;
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  const double *z = is3D() ? mZ.constData() : nullptr;
  const double *m = isMeasure() ? mM.constData() : nullptr;
 
  double prevX = *x++;
  double prevY = *y++;
  double prevZ = z ? *z++ : 0.0;
  double prevM = m ? *m++ : 0.0;
 
  if ( qgsDoubleNear( distance, 0.0 ) )
  {
    visitPoint( prevX, prevY, prevZ, prevM, prevX, prevY, prevZ, prevM, prevX, prevY, prevZ, prevM );
    return;
  }
 
  double pZ = std::numeric_limits<double>::quiet_NaN();
  double pM = std::numeric_limits<double>::quiet_NaN();
  double nextPointDistance = distance;
  const double eps = 4 * nextPointDistance * std::numeric_limits<double>::epsilon();
  for ( int i = 1; i < totalPoints; ++i )
  {
    double thisX = *x++;
    double thisY = *y++;
    double thisZ = z ? *z++ : 0.0;
    double thisM = m ? *m++ : 0.0;
 
    const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
    while ( nextPointDistance < distanceTraversed + segmentLength || qgsDoubleNear( nextPointDistance, distanceTraversed + segmentLength, eps ) )
    {
      // point falls on this segment - truncate to segment length if qgsDoubleNear test was actually > segment length
      const double distanceToPoint = std::min( nextPointDistance - distanceTraversed, segmentLength );
      double pX, pY;
      QgsGeometryUtilsBase::
        pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, pX, pY, z ? &prevZ : nullptr, z ? &thisZ : nullptr, z ? &pZ : nullptr, m ? &prevM : nullptr, m ? &thisM : nullptr, m ? &pM : nullptr );
 
      if ( !visitPoint( pX, pY, pZ, pM, prevX, prevY, prevZ, prevM, thisX, thisY, thisZ, thisM ) )
        return;
 
      nextPointDistance += distance;
    }
 
    distanceTraversed += segmentLength;
    prevX = thisX;
    prevY = thisY;
    prevZ = thisZ;
    prevM = thisM;
  }
}
 
QgsPoint *QgsLineString::interpolatePoint( const double distance ) const
{
  if ( distance < 0 )
    return nullptr;
 
  Qgis::WkbType pointType = Qgis::WkbType::Point;
  if ( is3D() )
    pointType = Qgis::WkbType::PointZ;
  if ( isMeasure() )
    pointType = QgsWkbTypes::addM( pointType );
 
  std::unique_ptr< QgsPoint > res;
  visitPointsByRegularDistance( distance, [&]( double x, double y, double z, double m, double, double, double, double, double, double, double, double ) -> bool {
    res = std::make_unique< QgsPoint >( pointType, x, y, z, m );
    return false;
  } );
  return res.release();
}
 
bool QgsLineString::lineLocatePointByM( double m, double &x, double &y, double &z, double &distanceFromStart, bool use3DDistance ) const
{
  return lineLocatePointByMPrivate( m, x, y, z, distanceFromStart, use3DDistance, false );
}
 
bool QgsLineString::lineLocatePointByMPrivate( double m, double &x, double &y, double &z, double &distanceFromStart, bool use3DDistance, bool haveInterpolatedM ) const
{
  if ( !isMeasure() )
    return false;
 
  distanceFromStart = 0;
  const int totalPoints = numPoints();
  if ( totalPoints == 0 )
    return false;
 
  const double *xData = mX.constData();
  const double *yData = mY.constData();
  const double *mData = mM.constData();
 
  const double *zData = is3D() ? mZ.constData() : nullptr;
  use3DDistance &= static_cast< bool >( zData );
 
  double prevX = *xData++;
  double prevY = *yData++;
  double prevZ = zData ? *zData++ : 0;
  double prevM = *mData++;
 
  int i = 1;
  while ( i < totalPoints )
  {
    double thisX = *xData++;
    double thisY = *yData++;
    double thisZ = zData ? *zData++ : 0;
    double thisM = *mData++;
    const double segmentLength = use3DDistance ? QgsGeometryUtilsBase::distance3D( thisX, thisY, thisZ, prevX, prevY, prevZ ) : QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
 
    if ( std::isnan( thisM ) )
    {
      if ( haveInterpolatedM )
        return false;
 
      // if we hit a NaN m value, interpolate m to fill the blanks and then re-try
      std::unique_ptr< QgsLineString > interpolatedM( interpolateM( use3DDistance ) );
      return interpolatedM->lineLocatePointByMPrivate( m, x, y, z, distanceFromStart, use3DDistance, true );
    }
    else
    {
      // check if target m value falls within this segment's range
      if ( ( prevM < m && thisM > m ) || ( prevM > m && thisM < m ) || qgsDoubleNear( prevM, m ) || qgsDoubleNear( thisM, m ) )
      {
        // use centroid for constant value m segments
        if ( qgsDoubleNear( thisM, m ) && ( i < totalPoints - 1 ) && qgsDoubleNear( *mData, m ) )
        {
          distanceFromStart += segmentLength;
          // scan ahead till we find a vertex with a different m
          double totalLengthOfSegmentsWithConstantM = 0;
          for ( int j = 0; j < ( totalPoints - i ); ++j )
          {
            if ( !qgsDoubleNear( *( mData + j ), m ) )
              break;
 
            const double segmentLength = use3DDistance ? QgsGeometryUtilsBase::distance3D( *( xData + j - 1 ), *( yData + j - 1 ), *( zData + j - 1 ), *( xData + j ), *( yData + j ), *( zData + j ) )
                                                       : QgsGeometryUtilsBase::distance2D( *( xData + j - 1 ), *( yData + j - 1 ), *( xData + j ), *( yData + j ) );
            totalLengthOfSegmentsWithConstantM += segmentLength;
          }
 
          distanceFromStart += totalLengthOfSegmentsWithConstantM / 2;
          std::unique_ptr< QgsPoint> point( interpolatePoint( distanceFromStart ) );
          if ( !point )
            return false;
          x = point->x();
          y = point->y();
          z = point->z();
          return true;
        }
 
        const double delta = ( m - prevM ) / ( thisM - prevM );
 
        const double distanceToPoint = delta * segmentLength;
 
        QgsGeometryUtilsBase::pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToPoint, x, y );
        z = prevZ + ( thisZ - prevZ ) * delta;
        distanceFromStart += distanceToPoint;
        return true;
      }
    }
 
    distanceFromStart += segmentLength;
    prevX = thisX;
    prevY = thisY;
    prevZ = thisZ;
    prevM = thisM;
    ++i;
  }
  return false;
}
 
QgsLineString *QgsLineString::curveSubstring( double startDistance, double endDistance ) const
{
  if ( startDistance < 0 && endDistance < 0 )
    return createEmptyWithSameType();
 
  endDistance = std::max( startDistance, endDistance );
 
  const int totalPoints = numPoints();
  if ( totalPoints == 0 )
    return clone();
 
  QVector< QgsPoint > substringPoints;
  substringPoints.reserve( totalPoints );
 
  Qgis::WkbType pointType = Qgis::WkbType::Point;
  if ( is3D() )
    pointType = Qgis::WkbType::PointZ;
  if ( isMeasure() )
    pointType = QgsWkbTypes::addM( pointType );
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  const double *z = is3D() ? mZ.constData() : nullptr;
  const double *m = isMeasure() ? mM.constData() : nullptr;
 
  double distanceTraversed = 0;
  double prevX = *x++;
  double prevY = *y++;
  double prevZ = z ? *z++ : 0.0;
  double prevM = m ? *m++ : 0.0;
  bool foundStart = false;
 
  if ( startDistance < 0 )
    startDistance = 0;
 
  for ( int i = 1; i < totalPoints; ++i )
  {
    double thisX = *x++;
    double thisY = *y++;
    double thisZ = z ? *z++ : 0.0;
    double thisM = m ? *m++ : 0.0;
 
    const double segmentLength = QgsGeometryUtilsBase::distance2D( thisX, thisY, prevX, prevY );
 
    if ( distanceTraversed <= startDistance && startDistance < distanceTraversed + segmentLength )
    {
      // start point falls on this segment
      const double distanceToStart = startDistance - distanceTraversed;
      double startX, startY;
      double startZ = 0;
      double startM = 0;
      QgsGeometryUtilsBase::
        pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToStart, startX, startY, z ? &prevZ : nullptr, z ? &thisZ : nullptr, z ? &startZ : nullptr, m ? &prevM : nullptr, m ? &thisM : nullptr, m ? &startM : nullptr );
      substringPoints << QgsPoint( pointType, startX, startY, startZ, startM );
      foundStart = true;
    }
    if ( foundStart && ( distanceTraversed + segmentLength > endDistance ) )
    {
      // end point falls on this segment
      const double distanceToEnd = endDistance - distanceTraversed;
      double endX, endY;
      double endZ = 0;
      double endM = 0;
      QgsGeometryUtilsBase::
        pointOnLineWithDistance( prevX, prevY, thisX, thisY, distanceToEnd, endX, endY, z ? &prevZ : nullptr, z ? &thisZ : nullptr, z ? &endZ : nullptr, m ? &prevM : nullptr, m ? &thisM : nullptr, m ? &endM : nullptr );
      substringPoints << QgsPoint( pointType, endX, endY, endZ, endM );
    }
    else if ( foundStart )
    {
      substringPoints << QgsPoint( pointType, thisX, thisY, thisZ, thisM );
    }
 
    prevX = thisX;
    prevY = thisY;
    prevZ = thisZ;
    prevM = thisM;
    distanceTraversed += segmentLength;
    if ( distanceTraversed >= endDistance )
      break;
  }
 
  // start point is the last node
  if ( !foundStart && qgsDoubleNear( distanceTraversed, startDistance ) )
  {
    substringPoints << QgsPoint( pointType, prevX, prevY, prevZ, prevM ) << QgsPoint( pointType, prevX, prevY, prevZ, prevM );
  }
 
  return new QgsLineString( substringPoints );
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
void QgsLineString::draw( QPainter &p ) const
{
  p.drawPolyline( asQPolygonF() );
}
 
void QgsLineString::addToPainterPath( QPainterPath &path ) const
{
  int nPoints = numPoints();
  if ( nPoints < 1 )
  {
    return;
  }
 
  if ( path.isEmpty() || path.currentPosition() != QPointF( mX.at( 0 ), mY.at( 0 ) ) )
  {
    path.moveTo( mX.at( 0 ), mY.at( 0 ) );
  }
 
  for ( int i = 1; i < nPoints; ++i )
  {
    path.lineTo( mX.at( i ), mY.at( i ) );
  }
}
 
void QgsLineString::drawAsPolygon( QPainter &p ) const
{
  p.drawPolygon( asQPolygonF() );
}
 
QgsCompoundCurve *QgsLineString::toCurveType() const
{
  QgsCompoundCurve *compoundCurve = new QgsCompoundCurve();
  compoundCurve->addCurve( clone() );
  return compoundCurve;
}
 
void QgsLineString::extend( double startDistance, double endDistance )
{
  if ( mX.size() < 2 || mY.size() < 2 )
    return;
 
  const bool extendStart = startDistance > 0;
  const bool extendEnd = endDistance > 0;
 
  // start of line
  if ( extendStart )
  {
    const double currentLen = std::sqrt( std::pow( mX.at( 0 ) - mX.at( 1 ), 2 ) + std::pow( mY.at( 0 ) - mY.at( 1 ), 2 ) );
    const double newLen = currentLen + startDistance;
    mX[0] = mX.at( 1 ) + ( mX.at( 0 ) - mX.at( 1 ) ) / currentLen * newLen;
    mY[0] = mY.at( 1 ) + ( mY.at( 0 ) - mY.at( 1 ) ) / currentLen * newLen;
  }
  // end of line
  if ( extendEnd )
  {
    const int last = mX.size() - 1;
    const double currentLen = std::sqrt( std::pow( mX.at( last ) - mX.at( last - 1 ), 2 ) + std::pow( mY.at( last ) - mY.at( last - 1 ), 2 ) );
    const double newLen = currentLen + endDistance;
    mX[last] = mX.at( last - 1 ) + ( mX.at( last ) - mX.at( last - 1 ) ) / currentLen * newLen;
    mY[last] = mY.at( last - 1 ) + ( mY.at( last ) - mY.at( last - 1 ) ) / currentLen * newLen;
  }
 
  if ( extendStart || extendEnd )
    clearCache(); //set bounding box invalid
}
 
QgsLineString *QgsLineString::createEmptyWithSameType() const
{
  auto result = std::make_unique< QgsLineString >();
  result->mWkbType = mWkbType;
  return result.release();
}
 
QString QgsLineString::geometryType() const
{
  return u"LineString"_s;
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
bool QgsLineString::insertVertex( QgsVertexId position, const QgsPoint &vertex )
{
  if ( position.vertex < 0 || position.vertex > mX.size() )
  {
    return false;
  }
 
  if ( mWkbType == Qgis::WkbType::Unknown || mX.isEmpty() )
  {
    setZMTypeFromSubGeometry( &vertex, Qgis::WkbType::LineString );
  }
 
  mX.insert( position.vertex, vertex.x() );
  mY.insert( position.vertex, vertex.y() );
  if ( is3D() )
  {
    mZ.insert( position.vertex, vertex.z() );
  }
  if ( isMeasure() )
  {
    mM.insert( position.vertex, vertex.m() );
  }
  clearCache(); //set bounding box invalid
  return true;
}
 
bool QgsLineString::deleteVertex( QgsVertexId position )
{
  if ( position.vertex >= mX.size() || position.vertex < 0 )
  {
    return false;
  }
 
  mX.remove( position.vertex );
  mY.remove( position.vertex );
  if ( is3D() )
  {
    mZ.remove( position.vertex );
  }
  if ( isMeasure() )
  {
    mM.remove( position.vertex );
  }
 
  if ( numPoints() == 1 )
  {
    clear();
  }
 
  clearCache(); //set bounding box invalid
  return true;
}
 
bool QgsLineString::deleteVertices( const QSet<QgsVertexId> &positions )
{
  if ( positions.isEmpty() )
  {
    return false;
  }
 
  QList<QgsVertexId> vertices( positions.begin(), positions.end() );
 
  for ( QgsVertexId pos : positions )
  {
    if ( !hasVertex( pos ) )
    {
      return false;
    }
  }
 
  if ( numPoints() <= vertices.size() + 1 )
  {
    clear();
    return true;
  }
 
  std::sort( vertices.begin(), vertices.end(), []( const QgsVertexId &a, const QgsVertexId &b ) { return a.vertex > b.vertex; } );
 
  for ( QgsVertexId position : vertices )
  {
    mX.remove( position.vertex );
    mY.remove( position.vertex );
    if ( is3D() )
    {
      mZ.remove( position.vertex );
    }
    if ( isMeasure() )
    {
      mM.remove( position.vertex );
    }
  }
 
  clearCache(); //set bounding box invalid
  return true;
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
void QgsLineString::addVertex( const QgsPoint &pt )
{
  if ( mWkbType == Qgis::WkbType::Unknown || mX.isEmpty() )
  {
    setZMTypeFromSubGeometry( &pt, Qgis::WkbType::LineString );
  }
 
  mX.append( pt.x() );
  mY.append( pt.y() );
  if ( is3D() )
  {
    mZ.append( pt.z() );
  }
  if ( isMeasure() )
  {
    mM.append( pt.m() );
  }
  clearCache(); //set bounding box invalid
}
 
double QgsLineString::closestSegment( const QgsPoint &pt, QgsPoint &segmentPt, QgsVertexId &vertexAfter, int *leftOf, double epsilon ) const
{
  double sqrDist = std::numeric_limits<double>::max();
  double leftOfDist = std::numeric_limits<double>::max();
  int prevLeftOf = 0;
  double prevLeftOfX = 0.0;
  double prevLeftOfY = 0.0;
  double testDist = 0;
  double segmentPtX, segmentPtY;
 
  if ( leftOf )
    *leftOf = 0;
 
  const int size = mX.size();
  if ( size == 0 || size == 1 )
  {
    vertexAfter = QgsVertexId( 0, 0, 0 );
    return -1;
  }
 
  const double *xData = mX.constData();
  const double *yData = mY.constData();
  for ( int i = 1; i < size; ++i )
  {
    double prevX = xData[i - 1];
    double prevY = yData[i - 1];
    double currentX = xData[i];
    double currentY = yData[i];
    testDist = QgsGeometryUtilsBase::sqrDistToLine( pt.x(), pt.y(), prevX, prevY, currentX, currentY, segmentPtX, segmentPtY, epsilon );
    if ( testDist < sqrDist )
    {
      sqrDist = testDist;
      segmentPt.setX( segmentPtX );
      segmentPt.setY( segmentPtY );
      vertexAfter.part = 0;
      vertexAfter.ring = 0;
      vertexAfter.vertex = i;
    }
    if ( leftOf && qgsDoubleNear( testDist, sqrDist ) )
    {
      int left = QgsGeometryUtilsBase::leftOfLine( pt.x(), pt.y(), prevX, prevY, currentX, currentY );
      // if left equals 0, the test could not be performed (e.g. point in line with segment or on segment)
      // so don't set leftOf in this case, and hope that there's another segment that's the same distance
      // where we can perform the check
      if ( left != 0 )
      {
        if ( qgsDoubleNear( testDist, leftOfDist ) && left != prevLeftOf && prevLeftOf != 0 )
        {
          // we have two possible segments each with equal distance to point, but they disagree
          // on whether or not the point is to the left of them.
          // so we test the segments themselves and flip the result.
          // see https://stackoverflow.com/questions/10583212/elegant-left-of-test-for-polyline
          *leftOf = -QgsGeometryUtilsBase::leftOfLine( currentX, currentY, prevLeftOfX, prevLeftOfY, prevX, prevY );
        }
        else
        {
          *leftOf = left;
        }
        prevLeftOf = *leftOf;
        leftOfDist = testDist;
        prevLeftOfX = prevX;
        prevLeftOfY = prevY;
      }
      else if ( testDist < leftOfDist )
      {
        *leftOf = left;
        leftOfDist = testDist;
        prevLeftOf = 0;
      }
    }
  }
  return sqrDist;
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
bool QgsLineString::pointAt( int node, QgsPoint &point, Qgis::VertexType &type ) const
{
  if ( node < 0 || node >= numPoints() )
  {
    return false;
  }
  point = pointN( node );
  type = Qgis::VertexType::Segment;
  return true;
}
 
QgsPoint QgsLineString::centroid() const
{
  if ( mX.isEmpty() )
    return QgsPoint();
 
  int numPoints = mX.count();
  if ( numPoints == 1 )
    return QgsPoint( mX.at( 0 ), mY.at( 0 ) );
 
  double totalLineLength = 0.0;
  double prevX = mX.at( 0 );
  double prevY = mY.at( 0 );
  double sumX = 0.0;
  double sumY = 0.0;
 
  for ( int i = 1; i < numPoints; ++i )
  {
    double currentX = mX.at( i );
    double currentY = mY.at( i );
    double segmentLength = std::sqrt( std::pow( currentX - prevX, 2.0 ) + std::pow( currentY - prevY, 2.0 ) );
    if ( qgsDoubleNear( segmentLength, 0.0 ) )
      continue;
 
    totalLineLength += segmentLength;
    sumX += segmentLength * ( currentX + prevX );
    sumY += segmentLength * ( currentY + prevY );
    prevX = currentX;
    prevY = currentY;
  }
  sumX *= 0.5;
  sumY *= 0.5;
 
  if ( qgsDoubleNear( totalLineLength, 0.0 ) )
    return QgsPoint( mX.at( 0 ), mY.at( 0 ) );
  else
    return QgsPoint( sumX / totalLineLength, sumY / totalLineLength );
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
void QgsLineString::sumUpArea( double &sum ) const
{
  if ( mHasCachedSummedUpArea )
  {
    sum += mSummedUpArea;
    return;
  }
 
  mSummedUpArea = 0;
  const int maxIndex = mX.size();
  if ( maxIndex < 2 )
  {
    mHasCachedSummedUpArea = true;
    return;
  }
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  double prevX = *x++;
  double prevY = *y++;
  for ( int i = 1; i < maxIndex; ++i )
  {
    mSummedUpArea += prevX * ( *y - prevY ) - prevY * ( *x - prevX );
    prevX = *x++;
    prevY = *y++;
  }
  mSummedUpArea *= 0.5;
 
  mHasCachedSummedUpArea = true;
  sum += mSummedUpArea;
}
 
void QgsLineString::sumUpArea3D( double &sum ) const
{
  if ( mHasCachedSummedUpArea3D )
  {
    sum += mSummedUpArea3D;
    return;
  }
 
  // No Z component. Fallback to the 2D version
  if ( mZ.isEmpty() )
  {
    double area2D = 0;
    sumUpArea( area2D );
    mSummedUpArea3D = area2D;
    mHasCachedSummedUpArea3D = true;
    sum += mSummedUpArea3D;
    return;
  }
 
  mSummedUpArea3D = 0;
 
  // Look for a reference unit normal
  QgsPoint ptA;
  QgsPoint ptB;
  QgsPoint ptC;
  if ( !QgsGeometryUtils::checkWeaklyFor3DPlane( this, ptA, ptB, ptC ) )
  {
    mHasCachedSummedUpArea3D = true;
    return;
  }
 
  QgsVector3D vAB = QgsVector3D( ptB.x() - ptA.x(), ptB.y() - ptA.y(), ptB.z() - ptA.z() );
  QgsVector3D vAC = QgsVector3D( ptC.x() - ptA.x(), ptC.y() - ptA.y(), ptC.z() - ptA.z() );
  QgsVector3D planeNormal = QgsVector3D::crossProduct( vAB, vAC );
 
  // Ensure a Consistent orientation: prioritize Z+, then Y+, then X+
  if ( !qgsDoubleNear( planeNormal.z(), 0.0 ) )
  {
    if ( planeNormal.z() < 0 )
    {
      planeNormal = -planeNormal;
    }
  }
  else if ( !qgsDoubleNear( planeNormal.y(), 0.0 ) )
  {
    if ( planeNormal.y() < 0 )
      planeNormal = -planeNormal;
  }
  else
  {
    if ( planeNormal.x() < 0 )
      planeNormal = -planeNormal;
  }
  planeNormal.normalize();
 
  const double *x = mX.constData();
  const double *y = mY.constData();
  const double *z = mZ.constData();
 
  double prevX = *x++;
  double prevY = *y++;
  double prevZ = *z++;
 
  double normalX = 0.0;
  double normalY = 0.0; // #spellok - Y component of normal vector
  double normalZ = 0.0;
 
  for ( unsigned int i = 1; i < mX.size(); ++i )
  {
    normalX += prevY * ( *z - prevZ ) - prevZ * ( *y - prevY );
    normalY += prevZ * ( *x - prevX ) - prevX * ( *z - prevZ ); // #spellok
    normalZ += prevX * ( *y - prevY ) - prevY * ( *x - prevX );
 
    prevX = *x++;
    prevY = *y++;
    prevZ = *z++;
  }
 
  mSummedUpArea3D = 0.5 * ( normalX * planeNormal.x() + normalY * planeNormal.y() + normalZ * planeNormal.z() ); // #spellok
 
  mHasCachedSummedUpArea3D = true;
  sum += mSummedUpArea3D;
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
void QgsLineString::close()
{
  if ( numPoints() < 1 || isClosed() )
  {
    return;
  }
  addVertex( startPoint() );
}
 
double QgsLineString::vertexAngle( QgsVertexId vertex ) const
{
  if ( mX.count() < 2 )
  {
    //undefined
    return 0.0;
  }
 
  if ( vertex.vertex == 0 || vertex.vertex >= ( numPoints() - 1 ) )
  {
    if ( isClosed() )
    {
      double previousX = mX.at( numPoints() - 2 );
      double previousY = mY.at( numPoints() - 2 );
      double currentX = mX.at( 0 );
      double currentY = mY.at( 0 );
      double afterX = mX.at( 1 );
      double afterY = mY.at( 1 );
      return QgsGeometryUtilsBase::averageAngle( previousX, previousY, currentX, currentY, afterX, afterY );
    }
    else if ( vertex.vertex == 0 )
    {
      return QgsGeometryUtilsBase::lineAngle( mX.at( 0 ), mY.at( 0 ), mX.at( 1 ), mY.at( 1 ) );
    }
    else
    {
      int a = numPoints() - 2;
      int b = numPoints() - 1;
      return QgsGeometryUtilsBase::lineAngle( mX.at( a ), mY.at( a ), mX.at( b ), mY.at( b ) );
    }
  }
  else
  {
    double previousX = mX.at( vertex.vertex - 1 );
    double previousY = mY.at( vertex.vertex - 1 );
    double currentX = mX.at( vertex.vertex );
    double currentY = mY.at( vertex.vertex );
    double afterX = mX.at( vertex.vertex + 1 );
    double afterY = mY.at( vertex.vertex + 1 );
    return QgsGeometryUtilsBase::averageAngle( previousX, previousY, currentX, currentY, afterX, afterY );
  }
}
 
double QgsLineString::segmentLength( QgsVertexId startVertex ) const
{
  if ( startVertex.vertex < 0 || startVertex.vertex >= mX.count() - 1 )
    return 0.0;
 
  double dx = mX.at( startVertex.vertex + 1 ) - mX.at( startVertex.vertex );
  double dy = mY.at( startVertex.vertex + 1 ) - mY.at( startVertex.vertex );
  return std::sqrt( dx * dx + dy * dy );
}
 
/***************************************************************************
 * This class is considered CRITICAL and any change MUST be accompanied with
 * full unit tests.
 * See details in QEP #17
 ****************************************************************************/
 
bool QgsLineString::convertTo( Qgis::WkbType type )
{
  if ( type == mWkbType )
    return true;
 
  clearCache();
  if ( type == Qgis::WkbType::LineString25D )
  {
    //special handling required for conversion to LineString25D
    dropMValue();
    addZValue( std::numeric_limits<double>::quiet_NaN() );
    mWkbType = Qgis::WkbType::LineString25D;
    return true;
  }
  else
  {
    return QgsCurve::convertTo( type );
  }
}
 
std::unique_ptr< QgsLineString > QgsLineString::measuredLine( double start, double end ) const
{
  const int nbpoints = numPoints();
  std::unique_ptr< QgsLineString > cloned( clone() );
 
  if ( !cloned->convertTo( QgsWkbTypes::addM( mWkbType ) ) )
  {
    return cloned;
  }
 
  if ( isEmpty() || ( nbpoints < 2 ) )
  {
    return cloned;
  }
 
  const double range = end - start;
  double lineLength = length();
  double lengthSoFar = 0.0;
 
 
  double *mOut = cloned->mM.data();
  *mOut++ = start;
  for ( int i = 1; i < nbpoints; ++i )
  {
    lengthSoFar += QgsGeometryUtilsBase::distance2D( mX[i - 1], mY[i - 1], mX[i], mY[i] );
    if ( lineLength > 0.0 )
      *mOut++ = start + range * lengthSoFar / lineLength;
    else if ( lineLength == 0.0 && nbpoints > 1 )
      *mOut++ = start + range * i / ( nbpoints - 1 );
    else
      *mOut++ = 0.0;
  }
 
  return cloned;
}
 
std::unique_ptr< QgsLineString > QgsLineString::interpolateM( bool use3DDistance ) const
{
  if ( !isMeasure() )
    return nullptr;
 
  const int totalPoints = numPoints();
  if ( totalPoints < 2 )
    return std::unique_ptr< QgsLineString >( clone() );
 
  const double *xData = mX.constData();
  const double *yData = mY.constData();
  const double *mData = mM.constData();
  const double *zData = is3D() ? mZ.constData() : nullptr;
  use3DDistance &= static_cast< bool >( zData );
 
  QVector< double > xOut( totalPoints );
  QVector< double > yOut( totalPoints );
  QVector< double > mOut( totalPoints );
  QVector< double > zOut( static_cast< bool >( zData ) ? totalPoints : 0 );
 
  double *xOutData = xOut.data();
  double *yOutData = yOut.data();
  double *mOutData = mOut.data();
  double *zOutData = static_cast< bool >( zData ) ? zOut.data() : nullptr;
 
  int i = 0;
  double currentSegmentLength = 0;
  double lastValidM = std::numeric_limits< double >::quiet_NaN();
  double prevX = *xData;
  double prevY = *yData;
  double prevZ = zData ? *zData : 0;
  while ( i < totalPoints )
  {
    double thisX = *xData++;
    double thisY = *yData++;
    double thisZ = zData ? *zData++ : 0;
    double thisM = *mData++;
 
    currentSegmentLength = use3DDistance ? QgsGeometryUtilsBase::distance3D( prevX, prevY, prevZ, thisX, thisY, thisZ ) : QgsGeometryUtilsBase::distance2D( prevX, prevY, thisX, thisY );
 
    if ( !std::isnan( thisM ) )
    {
      *xOutData++ = thisX;
      *yOutData++ = thisY;
      *mOutData++ = thisM;
      if ( zOutData )
        *zOutData++ = thisZ;
      lastValidM = thisM;
    }
    else if ( i == 0 )
    {
      // nan m value at start of line, read ahead to find first non-nan value and backfill
      int j = 0;
      double scanAheadM = thisM;
      while ( i + j + 1 < totalPoints && std::isnan( scanAheadM ) )
      {
        scanAheadM = mData[j];
        ++j;
      }
      if ( std::isnan( scanAheadM ) )
      {
        // no valid m values in line
        return nullptr;
      }
      *xOutData++ = thisX;
      *yOutData++ = thisY;
      *mOutData++ = scanAheadM;
      if ( zOutData )
        *zOutData++ = thisZ;
      for ( ; i < j; ++i )
      {
        thisX = *xData++;
        thisY = *yData++;
        *xOutData++ = thisX;
        *yOutData++ = thisY;
        *mOutData++ = scanAheadM;
        mData++;
        if ( zOutData )
          *zOutData++ = *zData++;
      }
      lastValidM = scanAheadM;
    }
    else
    {
      // nan m value in middle of line, read ahead till next non-nan value and interpolate
      int j = 0;
      double scanAheadX = thisX;
      double scanAheadY = thisY;
      double scanAheadZ = thisZ;
      double distanceToNextValidM = currentSegmentLength;
      std::vector< double > scanAheadSegmentLengths;
      scanAheadSegmentLengths.emplace_back( currentSegmentLength );
 
      double nextValidM = std::numeric_limits< double >::quiet_NaN();
      while ( i + j < totalPoints - 1 )
      {
        double nextScanAheadX = xData[j];
        double nextScanAheadY = yData[j];
        double nextScanAheadZ = zData ? zData[j] : 0;
        double nextScanAheadM = mData[j];
        const double scanAheadSegmentLength = use3DDistance ? QgsGeometryUtilsBase::distance3D( scanAheadX, scanAheadY, scanAheadZ, nextScanAheadX, nextScanAheadY, nextScanAheadZ )
                                                            : QgsGeometryUtilsBase::distance2D( scanAheadX, scanAheadY, nextScanAheadX, nextScanAheadY );
        scanAheadSegmentLengths.emplace_back( scanAheadSegmentLength );
        distanceToNextValidM += scanAheadSegmentLength;
 
        if ( !std::isnan( nextScanAheadM ) )
        {
          nextValidM = nextScanAheadM;
          break;
        }
 
        scanAheadX = nextScanAheadX;
        scanAheadY = nextScanAheadY;
        scanAheadZ = nextScanAheadZ;
        ++j;
      }
 
      if ( std::isnan( nextValidM ) )
      {
        // no more valid m values, so just fill remainder of vertices with previous valid m value
        *xOutData++ = thisX;
        *yOutData++ = thisY;
        *mOutData++ = lastValidM;
        if ( zOutData )
          *zOutData++ = thisZ;
        ++i;
        for ( ; i < totalPoints; ++i )
        {
          *xOutData++ = *xData++;
          *yOutData++ = *yData++;
          *mOutData++ = lastValidM;
          if ( zOutData )
            *zOutData++ = *zData++;
        }
        break;
      }
      else
      {
        // interpolate along segments
        const double delta = ( nextValidM - lastValidM ) / distanceToNextValidM;
        *xOutData++ = thisX;
        *yOutData++ = thisY;
        *mOutData++ = lastValidM + delta * scanAheadSegmentLengths[0];
        double totalScanAheadLength = scanAheadSegmentLengths[0];
        if ( zOutData )
          *zOutData++ = thisZ;
        for ( int k = 1; k <= j; ++i, ++k )
        {
          thisX = *xData++;
          thisY = *yData++;
          *xOutData++ = thisX;
          *yOutData++ = thisY;
          totalScanAheadLength += scanAheadSegmentLengths[k];
          *mOutData++ = lastValidM + delta * totalScanAheadLength;
          mData++;
          if ( zOutData )
            *zOutData++ = *zData++;
        }
        lastValidM = nextValidM;
      }
    }
 
    prevX = thisX;
    prevY = thisY;
    prevZ = thisZ;
    ++i;
  }
  return std::make_unique< QgsLineString >( xOut, yOut, zOut, mOut );
}
 
double QgsLineString::distanceBetweenVertices( QgsVertexId fromVertex, QgsVertexId toVertex ) const
{
  // Convert QgsVertexId to simple vertex numbers for linestrings (single ring, single part)
  if ( fromVertex.part != 0 || fromVertex.ring != 0 || toVertex.part != 0 || toVertex.ring != 0 )
    return -1.0;
 
  const int fromVertexNumber = fromVertex.vertex;
  const int toVertexNumber = toVertex.vertex;
 
  // Ensure fromVertex < toVertex for simplicity
  if ( fromVertexNumber > toVertexNumber )
  {
    return distanceBetweenVertices( QgsVertexId( 0, 0, toVertexNumber ), QgsVertexId( 0, 0, fromVertexNumber ) );
  }
 
  const int nPoints = numPoints();
  if ( fromVertexNumber < 0 || fromVertexNumber >= nPoints || toVertexNumber < 0 || toVertexNumber >= nPoints )
    return -1.0;
 
  if ( fromVertexNumber == toVertexNumber )
    return 0.0;
 
  const bool is3DGeometry = is3D();
  const double *xData = mX.constData();
  const double *yData = mY.constData();
  const double *zData = is3DGeometry ? mZ.constData() : nullptr;
  double totalDistance = 0.0;
 
  // For linestring, just accumulate Euclidean distances between consecutive points
  for ( int i = fromVertexNumber; i < toVertexNumber; ++i )
  {
    double dx = xData[i + 1] - xData[i];
    double dy = yData[i + 1] - yData[i];
    double dz = 0.0;
 
    if ( is3DGeometry )
    {
      dz = zData[i + 1] - zData[i];
    }
 
    totalDistance += std::sqrt( dx * dx + dy * dy + dz * dz );
  }
 
  return totalDistance;
}