qgsmeshtracerenderer.cpp

Go to the documentation of this file.

/***************************************************************************
                         qgsmeshtracerenderer.cpp
                         -------------------------
    begin                : November 2019
    copyright            : (C) 2019 by Vincent Cloarec
    email                : vcloarec at gmail dot com
 ***************************************************************************/
 
/***************************************************************************
 *                                                                         *
 *   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 "qgsmeshtracerenderer.h"
#include "qgsmeshlayerrenderer.h"
#include "qgsrendercontext.h"
#include "qgslinestring.h"
#include "qgsmeshlayerinterpolator.h"
#include "qgsmeshlayerutils.h"
 
#include <QPointer>
 
 
#ifndef M_DEG2RAD
#define M_DEG2RAD 0.0174532925
#endif
 
 
QgsVector QgsMeshVectorValueInterpolator::vectorValue( const QgsPointXY &point ) const
{
  if ( mCacheFaceIndex != -1 && mCacheFaceIndex < mTriangularMesh.triangles().count() )
  {
    QgsVector res = interpolatedValuePrivate( mCacheFaceIndex, point );
    if ( isVectorValid( res ) )
    {
      activeFaceFilter( res, mCacheFaceIndex );
      return res;
    }
  }
 
  //point is not on the face associated with mCacheIndex --> search for the face containing the point
  QList<int> potentialFaceIndexes = mTriangularMesh.faceIndexesForRectangle( QgsRectangle( point, point ) );
  mCacheFaceIndex = -1;
  for ( const int faceIndex : potentialFaceIndexes )
  {
    QgsVector res = interpolatedValuePrivate( faceIndex, point );
    if ( isVectorValid( res ) )
    {
      mCacheFaceIndex = faceIndex;
      activeFaceFilter( res, mCacheFaceIndex );
      return res;
    }
  }
 
  //--> no face found return non valid vector
  return ( QgsVector( std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN() ) );
 
}
 
QgsMeshVectorValueInterpolator &QgsMeshVectorValueInterpolator::operator=( const QgsMeshVectorValueInterpolator &other )
{
  mTriangularMesh = other.mTriangularMesh;
  mDatasetValues = other.mDatasetValues;
  mActiveFaceFlagValues = other.mActiveFaceFlagValues;
  mFaceCache = other.mFaceCache;
  mCacheFaceIndex = other.mCacheFaceIndex;
  mUseScalarActiveFaceFlagValues = other.mUseScalarActiveFaceFlagValues;
 
  return *this;
}
 
QgsMeshVectorValueInterpolatorFromVertex::
QgsMeshVectorValueInterpolatorFromVertex( const QgsTriangularMesh &triangularMesh, const QgsMeshDataBlock &datasetVectorValues )
  : QgsMeshVectorValueInterpolator( triangularMesh, datasetVectorValues )
{
 
}
 
QgsMeshVectorValueInterpolatorFromVertex::
QgsMeshVectorValueInterpolatorFromVertex( const QgsTriangularMesh &triangularMesh,
    const QgsMeshDataBlock &datasetVectorValues,
    const QgsMeshDataBlock &scalarActiveFaceFlagValues )
  : QgsMeshVectorValueInterpolator( triangularMesh, datasetVectorValues, scalarActiveFaceFlagValues )
{
 
}
 
QgsMeshVectorValueInterpolatorFromVertex::QgsMeshVectorValueInterpolatorFromVertex( const QgsMeshVectorValueInterpolatorFromVertex &other ):
  QgsMeshVectorValueInterpolator( other )
{}
 
QgsMeshVectorValueInterpolatorFromVertex *QgsMeshVectorValueInterpolatorFromVertex::clone()
{
  return new QgsMeshVectorValueInterpolatorFromVertex( *this );
}
 
QgsMeshVectorValueInterpolatorFromVertex &QgsMeshVectorValueInterpolatorFromVertex::
operator=( const QgsMeshVectorValueInterpolatorFromVertex &other )
{
  QgsMeshVectorValueInterpolator::operator=( other );
  return ( *this );
}
 
QgsVector QgsMeshVectorValueInterpolatorFromVertex::interpolatedValuePrivate( int faceIndex, const QgsPointXY point ) const
{
  QgsMeshFace face = mTriangularMesh.triangles().at( faceIndex );
 
  QgsPoint p1 = mTriangularMesh.vertices().at( face.at( 0 ) );
  QgsPoint p2 = mTriangularMesh.vertices().at( face.at( 1 ) );
  QgsPoint p3 = mTriangularMesh.vertices().at( face.at( 2 ) );
 
  QgsVector v1 = QgsVector( mDatasetValues.value( face.at( 0 ) ).x(),
                            mDatasetValues.value( face.at( 0 ) ).y() );
 
  QgsVector v2 = QgsVector( mDatasetValues.value( face.at( 1 ) ).x(),
                            mDatasetValues.value( face.at( 1 ) ).y() );
 
  QgsVector v3 = QgsVector( mDatasetValues.value( face.at( 2 ) ).x(),
                            mDatasetValues.value( face.at( 2 ) ).y() );
 
  return QgsMeshLayerUtils::interpolateVectorFromVerticesData(
           p1,
           p2,
           p3,
           v1,
           v2,
           v3,
           point );
}
 
QgsMeshVectorValueInterpolator::QgsMeshVectorValueInterpolator( const QgsTriangularMesh &triangularMesh,
    const QgsMeshDataBlock &datasetVectorValues ):
  mTriangularMesh( triangularMesh ),
  mDatasetValues( datasetVectorValues ),
  mUseScalarActiveFaceFlagValues( false )
{}
 
QgsMeshVectorValueInterpolator::QgsMeshVectorValueInterpolator( const QgsTriangularMesh &triangularMesh,
    const QgsMeshDataBlock &datasetVectorValues,
    const QgsMeshDataBlock &scalarActiveFaceFlagValues ):
  mTriangularMesh( triangularMesh ),
  mDatasetValues( datasetVectorValues ),
  mActiveFaceFlagValues( scalarActiveFaceFlagValues ),
  mUseScalarActiveFaceFlagValues( true )
{}
 
QgsMeshVectorValueInterpolator::QgsMeshVectorValueInterpolator( const QgsMeshVectorValueInterpolator &other ):
  mTriangularMesh( other.mTriangularMesh ),
  mDatasetValues( other.mDatasetValues ),
  mActiveFaceFlagValues( other.mActiveFaceFlagValues ),
  mFaceCache( other.mFaceCache ),
  mCacheFaceIndex( other.mCacheFaceIndex ),
  mUseScalarActiveFaceFlagValues( other.mUseScalarActiveFaceFlagValues )
{}
 
void QgsMeshVectorValueInterpolator::updateCacheFaceIndex( const QgsPointXY &point ) const
{
  if ( ! QgsMeshUtils::isInTriangleFace( point, mFaceCache, mTriangularMesh.vertices() ) )
  {
    mCacheFaceIndex = mTriangularMesh.faceIndexForPoint_v2( point );
  }
}
 
bool QgsMeshVectorValueInterpolator::isVectorValid( const QgsVector &v ) const
{
  return !( std::isnan( v.x() ) || std::isnan( v.y() ) );
 
}
 
void QgsMeshVectorValueInterpolator::activeFaceFilter( QgsVector &vector, int faceIndex ) const
{
  if ( mUseScalarActiveFaceFlagValues && ! mActiveFaceFlagValues.active( mTriangularMesh.trianglesToNativeFaces()[faceIndex] ) )
    vector = QgsVector( std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN() ) ;
}
 
QSize QgsMeshStreamField::size() const
{
  return mFieldSize;
}
 
QPoint QgsMeshStreamField::topLeft() const
{
  return mFieldTopLeftInDeviceCoordinates;
}
 
int QgsMeshStreamField::resolution() const
{
  return mFieldResolution;
}
 
QgsPointXY QgsMeshStreamField::positionToMapCoordinates( const QPoint &pixelPosition, const QgsPointXY &positionInPixel )
{
  QgsPointXY mapPoint = mMapToFieldPixel.toMapCoordinates( pixelPosition );
  mapPoint = mapPoint + QgsVector( positionInPixel.x() * mMapToFieldPixel.mapUnitsPerPixel(),
                                   positionInPixel.y() * mMapToFieldPixel.mapUnitsPerPixel() );
  return mapPoint;
}
 
QgsMeshStreamField::QgsMeshStreamField(
  const QgsTriangularMesh &triangularMesh,
  const QgsMeshDataBlock &dataSetVectorValues,
  const QgsMeshDataBlock &scalarActiveFaceFlagValues,
  const QgsRectangle &layerExtent,
  double magnitudeMaximum,
  bool dataIsOnVertices,
  const QgsRenderContext &rendererContext,
  const QgsInterpolatedLineColor &vectorColoring,
  int resolution ):
  mFieldResolution( resolution ),
  mVectorColoring( vectorColoring ),
  mRenderContext( rendererContext ),
  mLayerExtent( layerExtent ),
  mMaximumMagnitude( magnitudeMaximum )
{
  if ( dataIsOnVertices )
  {
    if ( scalarActiveFaceFlagValues.isValid() )
      mVectorValueInterpolator.reset( new QgsMeshVectorValueInterpolatorFromVertex(
                                        triangularMesh,
                                        dataSetVectorValues,
                                        scalarActiveFaceFlagValues ) );
    else
      mVectorValueInterpolator.reset( new QgsMeshVectorValueInterpolatorFromVertex(
                                        triangularMesh,
                                        dataSetVectorValues ) );
  }
  else
  {
    if ( scalarActiveFaceFlagValues.isValid() )
      mVectorValueInterpolator.reset( new QgsMeshVectorValueInterpolatorFromFace(
                                        triangularMesh,
                                        dataSetVectorValues,
                                        scalarActiveFaceFlagValues ) );
    else
      mVectorValueInterpolator.reset( new QgsMeshVectorValueInterpolatorFromFace(
                                        triangularMesh,
                                        dataSetVectorValues ) );
  }
}
 
QgsMeshStreamField::QgsMeshStreamField( const QgsMeshStreamField &other ):
  mFieldSize( other.mFieldSize ),
  mFieldResolution( other.mFieldResolution ),
  mPen( other.mPen ),
  mTraceImage( other.mTraceImage ),
  mMapToFieldPixel( other.mMapToFieldPixel ),
  mVectorColoring( other.mVectorColoring ),
  mDirectionField( other.mDirectionField ),
  mRenderContext( other.mRenderContext ),
  mPixelFillingCount( other.mPixelFillingCount ),
  mMaxPixelFillingCount( other.mMaxPixelFillingCount ),
  mLayerExtent( other.mLayerExtent ),
  mMapExtent( other.mMapExtent ),
  mFieldTopLeftInDeviceCoordinates( other.mFieldTopLeftInDeviceCoordinates ),
  mValid( other.mValid ),
  mMaximumMagnitude( other.mMaximumMagnitude ),
  mPixelFillingDensity( other.mPixelFillingDensity ),
  mMinMagFilter( other.mMinMagFilter ),
  mMinimizeFieldSize( other.mMinimizeFieldSize )
{
  mPainter.reset( new QPainter( &mTraceImage ) );
  mVectorValueInterpolator =
    std::unique_ptr<QgsMeshVectorValueInterpolator>( other.mVectorValueInterpolator->clone() );
}
 
QgsMeshStreamField::~QgsMeshStreamField()
{
  if ( mPainter )
    mPainter->end();
}
 
void QgsMeshStreamField::updateSize( const QgsRenderContext &renderContext )
{
  mMapExtent = renderContext.mapExtent();
  const QgsMapToPixel &deviceMapToPixel = renderContext.mapToPixel();
  QgsRectangle layerExtent;
  try
  {
    QgsCoordinateTransform extentTransform = renderContext.coordinateTransform();
    extentTransform.setBallparkTransformsAreAppropriate( true );
    layerExtent = extentTransform.transformBoundingBox( mLayerExtent );
  }
  catch ( QgsCsException &cse )
  {
    Q_UNUSED( cse )
    //if the transform fails, consider the whole map
    layerExtent = mMapExtent;
  }
 
  QgsRectangle interestZoneExtent;
  if ( mMinimizeFieldSize )
    interestZoneExtent = layerExtent.intersect( mMapExtent );
  else
    interestZoneExtent = mMapExtent;
 
  if ( interestZoneExtent == QgsRectangle() )
  {
    mValid = false;
    mFieldSize = QSize();
    mFieldTopLeftInDeviceCoordinates = QPoint();
    initField();
    return;
  }
 
  QgsRectangle fieldInterestZoneInDeviceCoordinates = QgsMeshLayerUtils::boundingBoxToScreenRectangle( deviceMapToPixel, interestZoneExtent );
  mFieldTopLeftInDeviceCoordinates =
    QPoint( static_cast<int>( std::round( fieldInterestZoneInDeviceCoordinates.xMinimum() ) ),
            static_cast<int>( std::round( fieldInterestZoneInDeviceCoordinates.yMinimum() ) ) );
  int fieldWidthInDeviceCoordinate = int( fieldInterestZoneInDeviceCoordinates.width() );
  int fieldHeightInDeviceCoordinate = int ( fieldInterestZoneInDeviceCoordinates.height() );
 
  int fieldWidth = int( fieldWidthInDeviceCoordinate / mFieldResolution );
  int fieldHeight = int( fieldHeightInDeviceCoordinate / mFieldResolution );
 
  //increase the field size if this size is not adjusted to extent of zone of interest in device coordinates
  if ( fieldWidthInDeviceCoordinate % mFieldResolution > 0 )
    fieldWidth++;
  if ( fieldHeightInDeviceCoordinate % mFieldResolution > 0 )
    fieldHeight++;
 
  if ( fieldWidth == 0 || fieldHeight == 0 )
  {
    mFieldSize = QSize();
    mOutputExtent = QgsRectangle();
  }
  else
  {
    mFieldSize.setWidth( fieldWidth );
    mFieldSize.setHeight( fieldHeight );
    QgsPointXY pt1 = deviceMapToPixel.toMapCoordinates( mFieldTopLeftInDeviceCoordinates );
    QgsPointXY pt2 = deviceMapToPixel.toMapCoordinates( mFieldTopLeftInDeviceCoordinates + QPoint( fieldWidth, fieldHeight ) );
    QgsPointXY pt3 = deviceMapToPixel.toMapCoordinates( mFieldTopLeftInDeviceCoordinates + QPoint( 0, fieldHeight ) );
    QgsPointXY pt4 = deviceMapToPixel.toMapCoordinates( mFieldTopLeftInDeviceCoordinates + QPoint( fieldWidth, 0 ) );
 
    mOutputExtent = QgsRectangle( std::min( {pt1.x(), pt2.x(), pt3.x(), pt4.x()} ),
                                  std::min( {pt1.y(), pt2.y(), pt3.y(), pt4.y()} ),
                                  std::max( {pt1.x(), pt2.x(), pt3.x(), pt4.x()} ),
                                  std::max( {pt1.y(), pt2.y(), pt3.y(), pt4.y()} ),
                                  true );
  }
 
  double mapUnitPerFieldPixel;
  if ( interestZoneExtent.width() > 0 )
    mapUnitPerFieldPixel = deviceMapToPixel.mapUnitsPerPixel() * mFieldResolution * mFieldSize.width() /
                           ( fieldWidthInDeviceCoordinate / static_cast<double>( mFieldResolution ) ) ;
  else
    mapUnitPerFieldPixel = 1e-8;
 
  int fieldRightDevice = mFieldTopLeftInDeviceCoordinates.x() + mFieldSize.width() * mFieldResolution;
  int fieldBottomDevice = mFieldTopLeftInDeviceCoordinates.y() + mFieldSize.height() * mFieldResolution;
  QgsPointXY fieldRightBottomMap = deviceMapToPixel.toMapCoordinates( fieldRightDevice, fieldBottomDevice );
 
  int fieldTopDevice = mFieldTopLeftInDeviceCoordinates.x();
  int fieldLeftDevice = mFieldTopLeftInDeviceCoordinates.y();
  QgsPointXY fieldTopLeftMap = deviceMapToPixel.toMapCoordinates( fieldTopDevice, fieldLeftDevice );
 
  double xc = ( fieldRightBottomMap.x() + fieldTopLeftMap.x() ) / 2;
  double yc = ( fieldTopLeftMap.y() + fieldRightBottomMap.y() ) / 2;
 
  mMapToFieldPixel = QgsMapToPixel( mapUnitPerFieldPixel,
                                    xc,
                                    yc,
                                    fieldWidth,
                                    fieldHeight,
                                    deviceMapToPixel.mapRotation()
                                  );
 
  initField();
  mValid = true;
}
 
void QgsMeshStreamField::updateSize( const QgsRenderContext &renderContext, int resolution )
{
  if ( renderContext.mapExtent() == mMapExtent && resolution == mFieldResolution )
    return;
  mFieldResolution = resolution;
 
  updateSize( renderContext );
}
 
bool QgsMeshStreamField::isValid() const
{
  return mValid;
}
 
void QgsMeshStreamField::addTrace( QgsPointXY startPoint )
{
  addTrace( mMapToFieldPixel.transform( startPoint ).toQPointF().toPoint() );
}
 
 
void QgsMeshStreamField::addRandomTraces()
{
  if ( mMaximumMagnitude > 0 )
    while ( ( mPixelFillingCount < mMaxPixelFillingCount ) &&
            mRenderContext.feedback() && !mRenderContext.feedback()->isCanceled() )
      addRandomTrace();
}
 
void QgsMeshStreamField::addRandomTrace()
{
  if ( !mValid )
    return;
 
  int xRandom =  1 + std::rand() / int( ( RAND_MAX + 1u ) / uint( mFieldSize.width() ) )  ;
  int yRandom = 1 + std::rand() / int ( ( RAND_MAX + 1u ) / uint( mFieldSize.height() ) ) ;
  addTrace( QPoint( xRandom, yRandom ) );
}
 
void QgsMeshStreamField::addGriddedTraces( int dx, int dy )
{
  int i = 0 ;
  while ( i < mFieldSize.width() && mRenderContext.feedback() && !mRenderContext.feedback()->isCanceled() )
  {
    int j = 0 ;
    while ( j < mFieldSize.height() && mRenderContext.feedback() && !mRenderContext.feedback()->isCanceled() )
    {
      addTrace( QPoint( i, j ) );
      j += dy;
    }
    i += dx;
  }
}
 
void QgsMeshStreamField::addTracesOnMesh( const QgsTriangularMesh &mesh, const QgsRectangle &extent )
{
  QList<int> facesInExtent = mesh.faceIndexesForRectangle( extent );
  QSet<int> vertices;
  for ( auto f : std::as_const( facesInExtent ) )
  {
    auto face = mesh.triangles().at( f );
    for ( auto i : std::as_const( face ) )
      vertices.insert( i );
  }
 
  for ( auto i : std::as_const( vertices ) )
  {
    addTrace( mesh.vertices().at( i ) );
  }
}
 
void QgsMeshStreamField::addTrace( QPoint startPixel )
{
  //This is where each traces are constructed
  if ( !mPainter )
    return;
 
  if ( isTraceExists( startPixel ) || isTraceOutside( startPixel ) )
    return;
 
  if ( !mVectorValueInterpolator )
    return;
 
  if ( !( mMaximumMagnitude > 0 ) )
    return;
 
  mPainter->setPen( mPen );
 
  //position in the pixelField
  double x1 = 0;
  double y1 = 0;
 
  std::list<QPair<QPoint, FieldData>> chunkTrace;
 
  QPoint currentPixel = startPixel;
  QgsVector vector;
  FieldData data;
  data.time = 1;
 
  while ( true )
  {
    QgsPointXY mapPosition = positionToMapCoordinates( currentPixel, QgsPointXY( x1, y1 ) );
    vector = mVectorValueInterpolator->vectorValue( mapPosition ) ;
 
    if ( std::isnan( vector.x() ) || std::isnan( vector.y() ) )
    {
      mPixelFillingCount++;
      setChunkTrace( chunkTrace );
      break;
    }
 
    /* nondimensional value :  Vu=2 when the particle need dt=1 to go through a pixel with the mMagMax magnitude
     * The nondimensional size of the side of a pixel is 2
     */
    vector = vector.rotateBy( -mMapToFieldPixel.mapRotation() * M_DEG2RAD );
    QgsVector vu = vector / mMaximumMagnitude * 2;
    data.magnitude = vector.length();
 
    double Vx = vu.x();
    double Vy = vu.y();
    double Vu = data.magnitude / mMaximumMagnitude * 2; //nondimensional vector magnitude
 
    if ( qgsDoubleNear( Vu, 0 ) )
    {
      // no trace anymore
      addPixelToChunkTrace( currentPixel, data, chunkTrace );
      simplifyChunkTrace( chunkTrace );
      setChunkTrace( chunkTrace );
      break;
    }
 
    //calculates where the particle will be after dt=1,
    QgsPointXY  nextPosition = QgsPointXY( x1, y1 ) + vu;
    int incX = 0;
    int incY = 0;
    if ( nextPosition.x() > 1 )
      incX = +1;
    if ( nextPosition.x() < -1 )
      incX = -1;
    if ( nextPosition.y() > 1 )
      incY = +1;
    if ( nextPosition.y() < -1 )
      incY = -1;
 
    if ( incX != 0 || incY != 0 )
    {
      data.directionX = incX;
      data.directionY = -incY;
      //the particule leave the current pixel --> store pixels, calculates where the particle is and change the current pixel
      if ( chunkTrace.empty() )
      {
        storeInField( QPair<QPoint, FieldData>( currentPixel, data ) );
      }
      if ( addPixelToChunkTrace( currentPixel, data, chunkTrace ) )
      {
        setChunkTrace( chunkTrace );
        clearChunkTrace( chunkTrace );
      }
 
      data.time = 1;
      currentPixel += QPoint( incX, -incY );
      x1 = nextPosition.x() - 2 * incX;
      y1 = nextPosition.y() - 2 * incY;
    }
    else
    {
      double x2, y2;
      /*the particule still in the pixel --> "push" the position with the vector value to join a border
       * and calculate the time spent to go to this border
       */
      if ( qgsDoubleNear( Vy, 0 ) )
      {
        y2 = y1;
        if ( Vx > 0 )
          incX = +1;
        else
          incX = -1;
 
        x2 = incX ;
      }
      else if ( qgsDoubleNear( Vx, 0 ) )
      {
        x2 = x1;
        if ( Vy > 0 )
          incY = +1;
        else
          incY = -1;
 
        y2 = incY ;
      }
      else
      {
        if ( Vy > 0 )
          x2 = x1 + ( 1 -  y1 ) * Vx / fabs( Vy ) ;
        else
          x2 = x1 + ( 1 + y1 ) * Vx / fabs( Vy ) ;
        if ( Vx > 0 )
          y2 = y1 + ( 1 - x1 ) * Vy / fabs( Vx ) ;
        else
          y2 = y1 + ( 1 + x1 ) * Vy / fabs( Vx ) ;
 
        if ( x2 >= 1 )
          x2 = 1;
 
        if ( x2 <= -1 )
          x2 = -1;
 
        if ( y2 >= 1 )
          y2 = 1;
 
        if ( y2 <= -1 )
          y2 = -1;
 
      }
 
      //calculate distance
      double dx = x2 - x1;
      double dy = y2 - y1;
      double dl = sqrt( dx * dx + dy * dy );
 
      data.time += static_cast<float>( dl / Vu ) ; //adimensional time step : this the time needed to go to the border of the pixel
      if ( data.time > 10000 ) //Guard to prevent that the particle never leave the pixel
      {
        addPixelToChunkTrace( currentPixel, data, chunkTrace );
        setChunkTrace( chunkTrace );
        break;
      }
      x1 = x2;
      y1 = y2;
    }
 
    //test if the new current pixel is already defined, if yes no need to continue
    if ( isTraceExists( currentPixel ) )
    {
      //Set the pixel in the chunk before adding the current pixel because this pixel is already defined
      setChunkTrace( chunkTrace );
      addPixelToChunkTrace( currentPixel, data, chunkTrace );
      break;
    }
 
    if ( isTraceOutside( currentPixel ) )
    {
      setChunkTrace( chunkTrace );
      break;
    }
 
    if ( mRenderContext.feedback() && mRenderContext.feedback()->isCanceled() )
      break;
 
    if ( mRenderContext.renderingStopped() )
      break;
  }
 
  drawTrace( startPixel );
}
 
void QgsMeshStreamField::setResolution( int width )
{
  mFieldResolution = width;
}
 
QSize QgsMeshStreamField::imageSize() const
{
  return mFieldSize * mFieldResolution;
}
 
QPointF QgsMeshStreamField::fieldToDevice( const QPoint &pixel ) const
{
  QPointF p( pixel );
  p = mFieldResolution * p + QPointF( mFieldResolution - 1, mFieldResolution - 1 ) / 2;
  return p;
}
 
bool QgsMeshStreamField::addPixelToChunkTrace( QPoint &pixel,
    QgsMeshStreamField::FieldData &data,
    std::list<QPair<QPoint, QgsMeshStreamField::FieldData> > &chunkTrace )
{
  chunkTrace.emplace_back( pixel, data );
  if ( chunkTrace.size() == 3 )
  {
    simplifyChunkTrace( chunkTrace );
    return true;
  }
  return false;
}
 
void QgsMeshStreamlinesField::initField()
{
  mField = QVector<bool>( mFieldSize.width() * mFieldSize.height(), false );
  mDirectionField = QVector<unsigned char>( mFieldSize.width() * mFieldSize.height(), static_cast<unsigned char>( int( 0 ) ) );
  initImage();
}
 
void QgsMeshStreamlinesField::initImage()
{
  mTraceImage = QImage();
  switch ( mVectorColoring.coloringMethod() )
  {
    case QgsInterpolatedLineColor::ColorRamp:
    {
      QSize imgSize = mFieldSize * mFieldResolution;
      QgsRenderContext fieldContext = mRenderContext;
 
      fieldContext.setMapToPixel( mMapToFieldPixel );
      std::unique_ptr<QgsMeshLayerInterpolator> mScalarInterpolator(
        new QgsMeshLayerInterpolator(
          mTriangularMesh,
          mMagValues,
          mScalarActiveFaceFlagValues,
          mDataType,
          fieldContext,
          imgSize ) );
 
      QgsRasterShader *sh = new QgsRasterShader();
      sh->setRasterShaderFunction( new QgsColorRampShader( mVectorColoring.colorRampShader() ) ); // takes ownership of fcn
      QgsSingleBandPseudoColorRenderer renderer( mScalarInterpolator.get(), 0, sh );  // takes ownership of sh
      if ( imgSize.isValid() )
      {
        std::unique_ptr<QgsRasterBlock> bl( renderer.block( 0, mOutputExtent, imgSize.width(), imgSize.height(), mFeedBack ) );
        mTraceImage = bl->image();
      }
    }
    break;
    case QgsInterpolatedLineColor::SingleColor:
    {
      mTraceImage = QImage( mFieldSize * mFieldResolution, QImage::Format_ARGB32_Premultiplied );
      QColor col = mVectorColoring.singleColor();
      mTraceImage.fill( col );
    }
    break;
  }
 
  if ( !mTraceImage.isNull() )
  {
    mPainter.reset( new QPainter( &mTraceImage ) );
    mPainter->setRenderHint( QPainter::Antialiasing, true );
 
    mDrawingTraceImage = QImage( mTraceImage.size(), QImage::Format_ARGB32_Premultiplied );
    mDrawingTraceImage.fill( Qt::transparent );
    mDrawingTracePainter.reset( new QPainter( &mDrawingTraceImage ) );
    mDrawingTracePainter->setRenderHint( QPainter::Antialiasing, true );
  }
}
 
void QgsMeshStreamField::clearChunkTrace( std::list<QPair<QPoint, QgsMeshStreamField::FieldData> > &chunkTrace )
{
  auto one_before_end = std::prev( chunkTrace.end() );
  chunkTrace.erase( chunkTrace.begin(), one_before_end );
}
 
void QgsMeshStreamField::simplifyChunkTrace( std::list<QPair<QPoint, FieldData> > &chunkTrace )
{
  if ( chunkTrace.size() != 3 )
    return;
 
  auto ip3 = chunkTrace.begin();
  auto ip1 = ip3++;
  auto ip2 = ip3++;
 
  while ( ip3 != chunkTrace.end() && ip2 != chunkTrace.end() )
  {
    QPoint v1 = ( *ip1 ).first - ( *ip2 ).first;
    QPoint v2 = ( *ip2 ).first - ( *ip3 ).first;
    if ( v1.x()*v2.x() + v1.y()*v2.y() == 0 )
    {
      ( *ip1 ).second.time += ( ( *ip2 ).second.time ) / 2;
      ( *ip3 ).second.time += ( ( *ip2 ).second.time ) / 2;
      ( *ip1 ).second.directionX += ( *ip2 ).second.directionX;
      ( *ip1 ).second.directionY += ( *ip2 ).second.directionY;
      chunkTrace.erase( ip2 );
    }
    ip1 = ip3++;
    ip2 = ip3++;
  }
}
 
QgsMeshStreamlinesField::QgsMeshStreamlinesField( const QgsTriangularMesh &triangularMesh,
    const QgsMeshDataBlock &datasetVectorValues,
    const QgsMeshDataBlock &scalarActiveFaceFlagValues,
    const QgsRectangle &layerExtent,
    double magMax,
    bool dataIsOnVertices,
    QgsRenderContext &rendererContext,
    const QgsInterpolatedLineColor vectorColoring )
  : QgsMeshStreamField(
      triangularMesh,
      datasetVectorValues,
      scalarActiveFaceFlagValues,
      layerExtent,
      magMax,
      dataIsOnVertices,
      rendererContext,
      vectorColoring )
  , mMagValues( QgsMeshLayerUtils::calculateMagnitudes( datasetVectorValues ) )
{
}
 
QgsMeshStreamlinesField::QgsMeshStreamlinesField(
  const QgsTriangularMesh &triangularMesh,
  const QgsMeshDataBlock &datasetVectorValues,
  const QgsMeshDataBlock &scalarActiveFaceFlagValues,
  const QVector<double> &datasetMagValues,
  const QgsRectangle &layerExtent,
  QgsMeshLayerRendererFeedback *feedBack,
  double magMax,
  bool dataIsOnVertices,
  QgsRenderContext &rendererContext,
  const QgsInterpolatedLineColor vectorColoring )
  : QgsMeshStreamField(
      triangularMesh,
      datasetVectorValues,
      scalarActiveFaceFlagValues,
      layerExtent,
      magMax,
      dataIsOnVertices,
      rendererContext,
      vectorColoring )
  , mTriangularMesh( triangularMesh )
  , mMagValues( datasetMagValues )
  , mScalarActiveFaceFlagValues( scalarActiveFaceFlagValues )
  , mDataType( dataIsOnVertices ? QgsMeshDatasetGroupMetadata::DataOnVertices : QgsMeshDatasetGroupMetadata::DataOnFaces )
  , mFeedBack( feedBack )
{
}
 
void QgsMeshStreamlinesField::compose()
{
  if ( !mPainter )
    return;
  mPainter->setCompositionMode( QPainter::CompositionMode_DestinationIn );
  mPainter->drawImage( 0, 0, mDrawingTraceImage );
}
 
void QgsMeshStreamlinesField::storeInField( const QPair<QPoint, FieldData> pixelData )
{
  int i = pixelData.first.x();
  int j = pixelData.first.y();
  if ( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() )
  {
    mField[j * mFieldSize.width() + i] = true;
    int d = pixelData.second.directionX + 2 + ( pixelData.second.directionY + 1 ) * 3;
    mDirectionField[j * mFieldSize.width() + i] = static_cast<unsigned char>( d );
  }
}
 
void QgsMeshStreamField::setChunkTrace( std::list<QPair<QPoint, FieldData> > &chunkTrace )
{
  auto p = chunkTrace.begin();
  while ( p != chunkTrace.end() )
  {
    storeInField( ( *p ) );
    mPixelFillingCount++;
    ++p;
  }
}
 
void QgsMeshStreamlinesField::drawTrace( const QPoint &start ) const
{
  if ( !isTraceExists( start ) || isTraceOutside( start ) )
    return;
 
  if ( !mDrawingTracePainter )
    return;
 
  QPoint pt1 = start;
  QPoint curPt = pt1;
  int fieldWidth = mFieldSize.width();
  QSet<QgsPointXY> path;
  unsigned char dir = 0;
  unsigned char prevDir = mDirectionField.at( pt1.y() * fieldWidth  + pt1.x() );
 
  QVector<double> xPoly;
  QVector<double> yPoly;
  QPointF devicePt = fieldToDevice( pt1 );
  xPoly.append( devicePt.x() );
  yPoly.append( devicePt.y() );
 
  while ( isTraceExists( curPt ) && !isTraceOutside( curPt ) && !path.contains( curPt ) )
  {
    dir = mDirectionField.at( curPt.y() * fieldWidth  + curPt.x() );
    if ( dir == 5 ) //no direction, static pixel
      break;
 
    const QPoint curPtDir( ( dir - 1 ) % 3 - 1, ( dir - 1 ) / 3 - 1 );
    const QPoint pt2 = curPt + curPtDir;
 
    if ( dir != prevDir )
    {
      path.insert( curPt );
      devicePt = fieldToDevice( curPt );
      xPoly.append( devicePt.x() );
      yPoly.append( devicePt.y() );
      prevDir = dir;
    }
    curPt = pt2;
  }
 
  if ( ! isTraceExists( curPt ) || isTraceOutside( curPt ) )
  {
    // just add the last point
    devicePt = fieldToDevice( curPt - QPoint( ( dir - 1 ) % 3 - 1, ( dir - 1 ) / 3 - 1 ) );
    xPoly.append( devicePt.x() );
    yPoly.append( devicePt.y() );
  }
 
  QgsGeometry geom( new QgsLineString( xPoly, yPoly ) );
  geom = geom.simplify( 1.5 * mFieldResolution ).smooth( 1, 0.25, -1.0, 45 );
  QPen pen = mPen;
  pen.setColor( QColor( 0, 0, 0, 255 ) );
  mDrawingTracePainter->setPen( pen );
  mDrawingTracePainter->drawPolyline( geom.asQPolygonF() );
}
 
bool QgsMeshStreamlinesField::isTraceExists( const QPoint &pixel ) const
{
  int i = pixel.x();
  int j = pixel.y();
  if ( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() )
  {
    return mField[j * mFieldSize.width() + i];
  }
 
  return false;
}
 
bool QgsMeshStreamField::isTraceOutside( const QPoint &pixel ) const
{
  int i = pixel.x();
  int j = pixel.y();
 
  return !( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() );
}
 
void QgsMeshStreamField::setMinimizeFieldSize( bool minimizeFieldSize )
{
  mMinimizeFieldSize = minimizeFieldSize;
}
 
QgsMeshStreamField &QgsMeshStreamField::operator=( const QgsMeshStreamField &other )
{
  mFieldSize = other.mFieldSize ;
  mFieldResolution = other.mFieldResolution;
  mPen = other.mPen;
  mTraceImage = other.mTraceImage ;
  mMapToFieldPixel = other.mMapToFieldPixel ;
  mOutputExtent = other.mOutputExtent;
  mVectorColoring = other.mVectorColoring;
  mDirectionField = other.mDirectionField;
  mRenderContext = other.mRenderContext;
  mPixelFillingCount = other.mPixelFillingCount ;
  mMaxPixelFillingCount = other.mMaxPixelFillingCount ;
  mLayerExtent = other.mLayerExtent ;
  mMapExtent = other.mMapExtent;
  mFieldTopLeftInDeviceCoordinates = other.mFieldTopLeftInDeviceCoordinates ;
  mValid = other.mValid ;
  mMaximumMagnitude = other.mMaximumMagnitude ;
  mPixelFillingDensity = other.mPixelFillingDensity ;
  mMinMagFilter = other.mMinMagFilter ;
  mMaxMagFilter = other.mMaxMagFilter ;
  mMinimizeFieldSize = other.mMinimizeFieldSize ;
  mVectorValueInterpolator =
    std::unique_ptr<QgsMeshVectorValueInterpolator>( other.mVectorValueInterpolator->clone() );
 
  mPainter.reset( new QPainter( &mTraceImage ) );
 
  return ( *this );
}
 
void QgsMeshStreamField::initImage()
{
  mTraceImage = QImage( mFieldSize * mFieldResolution, QImage::Format_ARGB32 );
  if ( !mTraceImage.isNull() )
  {
    mTraceImage.fill( 0X00000000 );
    mPainter.reset( new QPainter( &mTraceImage ) );
    mPainter->setRenderHint( QPainter::Antialiasing, true );
    mPainter->setPen( mPen );
  }
}
 
bool QgsMeshStreamField::filterMag( double value ) const
{
  return ( mMinMagFilter < 0 || value > mMinMagFilter ) && ( mMaxMagFilter < 0 || value < mMaxMagFilter );
}
 
QImage QgsMeshStreamField::image() const
{
  if ( mTraceImage.isNull() )
    return QImage();
  return mTraceImage.scaled( mFieldSize * mFieldResolution, Qt::IgnoreAspectRatio, Qt::SmoothTransformation );
}
 
void QgsMeshStreamField::setPixelFillingDensity( double maxFilling )
{
  mPixelFillingDensity = maxFilling;
  mMaxPixelFillingCount = int( mPixelFillingDensity * mFieldSize.width() * mFieldSize.height() );
}
 
void QgsMeshStreamField::setColor( QColor color )
{
  mPen.setColor( color );
}
 
void QgsMeshStreamField::setLineWidth( double width )
{
  mPen.setWidthF( width );
}
 
void QgsMeshStreamField::setFilter( double min, double max )
{
  mMinMagFilter = min;
  mMaxMagFilter = max;
}
 
QgsMeshVectorValueInterpolatorFromFace::QgsMeshVectorValueInterpolatorFromFace( const QgsTriangularMesh &triangularMesh, const QgsMeshDataBlock &datasetVectorValues ):
  QgsMeshVectorValueInterpolator( triangularMesh, datasetVectorValues )
{}
 
QgsMeshVectorValueInterpolatorFromFace::QgsMeshVectorValueInterpolatorFromFace( const QgsTriangularMesh &triangularMesh, const QgsMeshDataBlock &datasetVectorValues, const QgsMeshDataBlock &scalarActiveFaceFlagValues ):
  QgsMeshVectorValueInterpolator( triangularMesh, datasetVectorValues, scalarActiveFaceFlagValues )
{}
 
QgsMeshVectorValueInterpolatorFromFace::QgsMeshVectorValueInterpolatorFromFace( const QgsMeshVectorValueInterpolatorFromFace &other ):
  QgsMeshVectorValueInterpolator( other )
{}
 
QgsMeshVectorValueInterpolatorFromFace *QgsMeshVectorValueInterpolatorFromFace::clone()
{
  return new QgsMeshVectorValueInterpolatorFromFace( *this );
}
 
QgsMeshVectorValueInterpolatorFromFace &QgsMeshVectorValueInterpolatorFromFace::operator=( const QgsMeshVectorValueInterpolatorFromFace &other )
{
  QgsMeshVectorValueInterpolator::operator=( other );
  return ( *this );
}
 
QgsVector QgsMeshVectorValueInterpolatorFromFace::interpolatedValuePrivate( int faceIndex, const QgsPointXY point ) const
{
  QgsMeshFace face = mTriangularMesh.triangles().at( faceIndex );
 
  QgsPoint p1 = mTriangularMesh.vertices().at( face.at( 0 ) );
  QgsPoint p2 = mTriangularMesh.vertices().at( face.at( 1 ) );
  QgsPoint p3 = mTriangularMesh.vertices().at( face.at( 2 ) );
 
  QgsVector vect = QgsVector( mDatasetValues.value( mTriangularMesh.trianglesToNativeFaces().at( faceIndex ) ).x(),
                              mDatasetValues.value( mTriangularMesh.trianglesToNativeFaces().at( faceIndex ) ).y() );
 
  return QgsMeshLayerUtils::interpolateVectorFromFacesData(
           p1,
           p2,
           p3,
           vect,
           point );
}
 
QgsMeshVectorStreamlineRenderer::QgsMeshVectorStreamlineRenderer(
  const QgsTriangularMesh &triangularMesh,
  const QgsMeshDataBlock &dataSetVectorValues,
  const QgsMeshDataBlock &scalarActiveFaceFlagValues,
  bool dataIsOnVertices,
  const QgsMeshRendererVectorSettings &settings,
  QgsRenderContext &rendererContext,
  const QgsRectangle &layerExtent, double magMax ):
  QgsMeshVectorStreamlineRenderer(
    triangularMesh,
    dataSetVectorValues,
    scalarActiveFaceFlagValues,
    QgsMeshLayerUtils::calculateMagnitudes( dataSetVectorValues ),
    dataIsOnVertices,
    settings, rendererContext,
    layerExtent,
    nullptr,
    magMax )
{}
 
QgsMeshVectorStreamlineRenderer::QgsMeshVectorStreamlineRenderer( const QgsTriangularMesh &triangularMesh,
    const QgsMeshDataBlock &dataSetVectorValues,
    const QgsMeshDataBlock &scalarActiveFaceFlagValues,
    const QVector<double> &datasetMagValues,
    bool dataIsOnVertices,
    const QgsMeshRendererVectorSettings &settings,
    QgsRenderContext &rendererContext,
    const QgsRectangle &layerExtent, QgsMeshLayerRendererFeedback *feedBack,
    double magMax ):
  mRendererContext( rendererContext )
{
  mStreamlineField.reset(
    new QgsMeshStreamlinesField(
      triangularMesh,
      dataSetVectorValues,
      scalarActiveFaceFlagValues,
      datasetMagValues,
      layerExtent,
      feedBack,
      magMax,
      dataIsOnVertices,
      rendererContext,
      settings.vectorStrokeColoring() ) );
 
  mStreamlineField->updateSize( rendererContext );
  mStreamlineField->setPixelFillingDensity( settings.streamLinesSettings().seedingDensity() );
  mStreamlineField->setLineWidth( rendererContext.convertToPainterUnits( settings.lineWidth(),
                                  Qgis::RenderUnit::Millimeters ) ) ;
  mStreamlineField->setColor( settings.color() );
 
  mStreamlineField->setFilter( settings.filterMin(), settings.filterMax() );
 
  switch ( settings.streamLinesSettings().seedingMethod() )
  {
    case QgsMeshRendererVectorStreamlineSettings::MeshGridded:
      if ( settings.isOnUserDefinedGrid() )
        mStreamlineField->addGriddedTraces( settings.userGridCellWidth(), settings.userGridCellHeight() );
      else
        mStreamlineField->addTracesOnMesh( triangularMesh, rendererContext.mapExtent() );
      break;
    case QgsMeshRendererVectorStreamlineSettings::Random:
      mStreamlineField->addRandomTraces();
      break;
  }
}
 
void QgsMeshVectorStreamlineRenderer::draw()
{
  if ( mRendererContext.renderingStopped() )
    return;
  mStreamlineField->compose();
  mRendererContext.painter()->drawImage( mStreamlineField->topLeft(), mStreamlineField->image() );
}
 
QgsMeshParticleTracesField::QgsMeshParticleTracesField( const QgsTriangularMesh &triangularMesh,
    const QgsMeshDataBlock &datasetVectorValues,
    const QgsMeshDataBlock &scalarActiveFaceFlagValues,
    const QgsRectangle &layerExtent,
    double magMax,
    bool dataIsOnVertices,
    const QgsRenderContext &rendererContext,
    const QgsInterpolatedLineColor vectorColoring ):
  QgsMeshStreamField( triangularMesh,
                      datasetVectorValues,
                      scalarActiveFaceFlagValues,
                      layerExtent,
                      magMax,
                      dataIsOnVertices,
                      rendererContext,
                      vectorColoring )
{
  std::srand( uint( ::time( nullptr ) ) );
  mPen.setCapStyle( Qt::RoundCap );
}
 
QgsMeshParticleTracesField::QgsMeshParticleTracesField( const QgsMeshParticleTracesField &other )
  : QgsMeshStreamField( other )
  , mTimeField( other.mTimeField )
  , mMagnitudeField( other.mMagnitudeField )
  , mParticles( other.mParticles )
  , mStumpImage( other.mStumpImage )
  , mTimeStep( other.mTimeStep )
  , mParticlesLifeTime( other.mParticlesLifeTime )
  , mParticlesCount( other.mParticlesCount )
  , mTailFactor( other.mTailFactor )
  , mParticleColor( other.mParticleColor )
  , mParticleSize( other.mParticleSize )
  , mStumpFactor( other.mStumpFactor )
{}
 
void QgsMeshParticleTracesField::addParticle( const QPoint &startPoint, double lifeTime )
{
  addTrace( startPoint );
  if ( time( startPoint ) > 0 )
  {
    QgsMeshTraceParticle p;
    p.lifeTime = lifeTime;
    p.position = startPoint;
    mParticles.append( p );
  }
 
}
 
void QgsMeshParticleTracesField::addParticleXY( const QgsPointXY &startPoint, double lifeTime )
{
  addParticle( mMapToFieldPixel.transform( startPoint ).toQPointF().toPoint(), lifeTime );
}
 
void QgsMeshParticleTracesField::moveParticles()
{
  stump();
  for ( auto &p : mParticles )
  {
    double spentTime = p.remainingTime; //adjust with the past remaining time
    size_t countAdded = 0;
    while ( spentTime < mTimeStep && p.lifeTime > 0 )
    {
      double timeToSpend = double( time( p.position ) );
      if ( timeToSpend > 0 )
      {
        p.lifeTime -= timeToSpend;
        spentTime += timeToSpend;
        QPoint dir = direction( p.position );
        if ( p.lifeTime > 0 )
        {
          p.position += dir;
          p.tail.emplace_back( p.position );
          countAdded++;
        }
        else
        {
          break;
        }
      }
      else
      {
        p.lifeTime = -1;
        break;
      }
    }
 
    if ( p.lifeTime <= 0 )
    {
      // the particle is not alive anymore
      p.lifeTime = 0;
      p.tail.clear();
    }
    else
    {
      p.remainingTime = spentTime - mTimeStep;
      while ( static_cast<int>( p.tail.size() )  >  mMinTailLength &&
              static_cast<double>( p.tail.size() ) > ( static_cast<double>( countAdded ) * mTailFactor ) )
        p.tail.erase( p.tail.begin() );
      drawParticleTrace( p );
    }
  }
 
  //remove empty (dead particles)
  int i = 0;
  while ( i < mParticles.count() )
  {
    if ( mParticles.at( i ).tail.size() == 0 )
      mParticles.removeAt( i );
    else
      ++i;
  }
 
  //add new particles if needed
  if ( mParticles.count() < mParticlesCount )
    addRandomParticles();
}
 
void QgsMeshParticleTracesField::addRandomParticles()
{
  if ( !isValid() )
    return;
 
  if ( mParticlesCount < 0 ) //for tests, add one particle on the center of the map
  {
    addParticleXY( QgsPointXY( mMapToFieldPixel.xCenter(), mMapToFieldPixel.yCenter() ), mParticlesLifeTime );
    return;
  }
 
  int count = mParticlesCount - mParticles.count();
 
  for ( int i = 0; i < count; ++i )
  {
    int xRandom =  1 + std::rand() / int( ( RAND_MAX + 1u ) / uint( mFieldSize.width() ) )  ;
    int yRandom = 1 + std::rand() / int ( ( RAND_MAX + 1u ) / uint( mFieldSize.height() ) ) ;
    double lifeTime = ( std::rand() / ( ( RAND_MAX + 1u ) / mParticlesLifeTime ) );
    addParticle( QPoint( xRandom, yRandom ), lifeTime );
  }
}
 
void QgsMeshParticleTracesField::storeInField( const QPair<QPoint, QgsMeshStreamField::FieldData> pixelData )
{
  int i = pixelData.first.x();
  int j = pixelData.first.y();
  if ( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() )
  {
    mTimeField[j * mFieldSize.width() + i] = pixelData.second.time;
    int d = pixelData.second.directionX + 2 + ( pixelData.second.directionY + 1 ) * 3;
    mDirectionField[j * mFieldSize.width() + i] = static_cast<unsigned char>( d );
    mMagnitudeField[j * mFieldSize.width() + i] = static_cast<float>( pixelData.second.magnitude );
  }
}
 
void QgsMeshParticleTracesField::initField()
{
  mTimeField = QVector<float>( mFieldSize.width() * mFieldSize.height(), -1 );
  mDirectionField = QVector<unsigned char>( mFieldSize.width() * mFieldSize.height(), static_cast<unsigned char>( int( 0 ) ) );
  mMagnitudeField = QVector<float>( mFieldSize.width() * mFieldSize.height(), 0 );
  initImage();
  mStumpImage = QImage( mFieldSize * mFieldResolution, QImage::Format_ARGB32 );
  mStumpImage.fill( QColor( 0, 0, 0, mStumpFactor ) ); //alpha=0 -> no persitence, alpha=255 -> total persistence
}
 
bool QgsMeshParticleTracesField::isTraceExists( const QPoint &pixel ) const
{
  int i = pixel.x();
  int j = pixel.y();
  if ( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() )
  {
    return mTimeField[j * mFieldSize.width() + i] >= 0;
  }
 
  return false;
}
 
void QgsMeshParticleTracesField::setStumpParticleWithLifeTime( bool stumpParticleWithLifeTime )
{
  mStumpParticleWithLifeTime = stumpParticleWithLifeTime;
}
 
void QgsMeshParticleTracesField::setParticlesColor( const QColor &c )
{
  mVectorColoring.setColor( c );
}
 
QgsMeshParticleTracesField &QgsMeshParticleTracesField::operator=( const QgsMeshParticleTracesField &other )
{
  QgsMeshStreamField::operator=( other );
  mTimeField = other.mTimeField;
  mDirectionField = other.mDirectionField;
  mParticles = other.mParticles;
  mStumpImage = other.mStumpImage;
  mTimeStep = other.mTimeStep;
  mParticlesLifeTime = other.mParticlesLifeTime;
  mParticlesCount = other.mParticlesCount;
  mMinTailLength = other.mMinTailLength;
  mTailFactor = other.mTailFactor;
  mParticleColor = other.mParticleColor;
  mParticleSize = other.mParticleSize;
  mStumpFactor = other.mStumpFactor;
  mStumpParticleWithLifeTime = other.mStumpParticleWithLifeTime;
 
  return ( *this );
}
 
void QgsMeshParticleTracesField::setMinTailLength( int minTailLength )
{
  mMinTailLength = minTailLength;
}
 
void QgsMeshParticleTracesField::setTailFactor( double tailFactor )
{
  mTailFactor = tailFactor;
}
 
void QgsMeshParticleTracesField::setParticleSize( double particleSize )
{
  mParticleSize = particleSize;
}
 
void QgsMeshParticleTracesField::setTimeStep( double timeStep )
{
  mTimeStep = timeStep;
}
 
void QgsMeshParticleTracesField::setParticlesLifeTime( double particlesLifeTime )
{
  mParticlesLifeTime = particlesLifeTime;
}
 
QImage QgsMeshParticleTracesField::imageRendered() const
{
  return mTraceImage;
}
 
void QgsMeshParticleTracesField::stump()
{
  if ( !mPainter )
    return;
  QgsScopedQPainterState painterState( mPainter.get() );
  mPainter->setCompositionMode( QPainter::CompositionMode_DestinationIn );
  mPainter->drawImage( QPoint( 0, 0 ), mStumpImage );
}
 
void QgsMeshParticleTracesField::setStumpFactor( int sf )
{
  mStumpFactor = sf;
  mStumpImage = QImage( mFieldSize * mFieldResolution, QImage::Format_ARGB32 );
  mStumpImage.fill( QColor( 0, 0, 0, mStumpFactor ) );
}
 
QPoint QgsMeshParticleTracesField::direction( QPoint position ) const
{
  int i = position.x();
  int j = position.y();
  if ( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() )
  {
    int dir = static_cast<int>( mDirectionField[j * mFieldSize.width() + i] );
    if ( dir != 0 && dir < 10 )
      return QPoint( ( dir - 1 ) % 3 - 1, ( dir - 1 ) / 3 - 1 );
  }
  return QPoint( 0, 0 );
}
 
float QgsMeshParticleTracesField::time( QPoint position ) const
{
  int i = position.x();
  int j = position.y();
  if ( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() )
  {
    return mTimeField[j * mFieldSize.width() + i];
  }
  return -1;
}
 
float QgsMeshParticleTracesField::magnitude( QPoint position ) const
{
  int i = position.x();
  int j = position.y();
  if ( i >= 0 && i < mFieldSize.width() && j >= 0 && j < mFieldSize.height() )
  {
    return mMagnitudeField[j * mFieldSize.width() + i];
  }
  return -1;
}
 
void QgsMeshParticleTracesField::drawParticleTrace( const QgsMeshTraceParticle &particle )
{
  if ( !mPainter )
    return;
  const std::list<QPoint> &tail = particle.tail;
  if ( tail.size() == 0 )
    return;
  double iniWidth = mParticleSize;
 
  size_t pixelCount = tail.size();
 
  double transparency = 1;
  if ( mStumpParticleWithLifeTime )
    transparency = sin( M_PI * particle.lifeTime / mParticlesLifeTime );
 
  double dw;
  if ( pixelCount > 1 )
    dw = iniWidth   / static_cast<double>( pixelCount );
  else
    dw = 0;
 
  auto ip1 = std::prev( tail.end() );
  auto ip2 = std::prev( ip1 );
  int i = 0;
  while ( ip1 != tail.begin() )
  {
    QPointF p1 = fieldToDevice( ( *ip1 ) );
    QPointF p2 = fieldToDevice( ( *ip2 ) );
    QColor traceColor = mVectorColoring.color( magnitude( *ip1 ) );
    traceColor.setAlphaF( traceColor.alphaF()*transparency );
    mPen.setColor( traceColor );
    mPen.setWidthF( iniWidth - i * dw );
    mPainter->setPen( mPen );
    mPainter->drawLine( p1, p2 );
    ip1--;
    ip2--;
    ++i;
  }
}
 
void QgsMeshParticleTracesField::setParticlesCount( int particlesCount )
{
  mParticlesCount = particlesCount;
}
 
QgsMeshVectorTraceAnimationGenerator::QgsMeshVectorTraceAnimationGenerator( const QgsTriangularMesh &triangularMesh,
    const QgsMeshDataBlock &dataSetVectorValues,
    const QgsMeshDataBlock &scalarActiveFaceFlagValues,
    bool dataIsOnVertices,
    const QgsRenderContext &rendererContext,
    const QgsRectangle &layerExtent,
    double magMax,
    const QgsMeshRendererVectorSettings &vectorSettings )
  : mParticleField( new QgsMeshParticleTracesField(
                      triangularMesh,
                      dataSetVectorValues,
                      scalarActiveFaceFlagValues,
                      layerExtent,
                      magMax,
                      dataIsOnVertices,
                      rendererContext,
                      vectorSettings.vectorStrokeColoring() ) )
  , mRendererContext( rendererContext )
{
  mParticleField->updateSize( rendererContext ) ;
}
 
QgsMeshVectorTraceAnimationGenerator::QgsMeshVectorTraceAnimationGenerator( QgsMeshLayer *layer, const QgsRenderContext &rendererContext ):
  mRendererContext( rendererContext )
{
  if ( !layer->triangularMesh() )
    layer->reload();
 
  QgsMeshDataBlock vectorDatasetValues;
  QgsMeshDataBlock scalarActiveFaceFlagValues;
  bool vectorDataOnVertices;
  double magMax;
 
  QgsMeshDatasetIndex datasetIndex = layer->activeVectorDatasetAtTime( rendererContext.temporalRange() );
 
  // Find out if we can use cache up to date. If yes, use it and return
  int datasetGroupCount = layer->dataProvider()->datasetGroupCount();
  const QgsMeshRendererVectorSettings vectorSettings = layer->rendererSettings().vectorSettings( datasetIndex.group() );
  QgsMeshLayerRendererCache *cache = layer->rendererCache();
 
  if ( ( cache->mDatasetGroupsCount == datasetGroupCount ) &&
       ( cache->mActiveVectorDatasetIndex == datasetIndex ) )
  {
    vectorDatasetValues = cache->mVectorDatasetValues;
    scalarActiveFaceFlagValues = cache->mScalarActiveFaceFlagValues;
    magMax = cache->mVectorDatasetMagMaximum;
    vectorDataOnVertices = cache->mVectorDataType == QgsMeshDatasetGroupMetadata::DataOnVertices;
  }
  else
  {
    const QgsMeshDatasetGroupMetadata metadata =
      layer->dataProvider()->datasetGroupMetadata( datasetIndex.group() );
    magMax = metadata.maximum();
    vectorDataOnVertices = metadata.dataType() == QgsMeshDatasetGroupMetadata::DataOnVertices;
 
    int count;
    if ( vectorDataOnVertices )
      count = layer->nativeMesh()->vertices.count();
    else
      count = layer->nativeMesh()->faces.count();
 
    vectorDatasetValues = QgsMeshLayerUtils::datasetValues( layer, datasetIndex, 0, count );
 
    scalarActiveFaceFlagValues = layer->dataProvider()->areFacesActive(
                                   datasetIndex,
                                   0,
                                   layer->nativeMesh()->faces.count() );
  }
 
  mParticleField = std::unique_ptr<QgsMeshParticleTracesField>( new QgsMeshParticleTracesField( ( *layer->triangularMesh() ),
                   vectorDatasetValues,
                   scalarActiveFaceFlagValues,
                   layer->extent(),
                   magMax,
                   vectorDataOnVertices,
                   rendererContext,
                   vectorSettings.vectorStrokeColoring() ) )  ;
 
  mParticleField->setMinimizeFieldSize( false );
  mParticleField->updateSize( mRendererContext );
}
 
QgsMeshVectorTraceAnimationGenerator::QgsMeshVectorTraceAnimationGenerator( const QgsMeshVectorTraceAnimationGenerator &other )
  : mParticleField( new QgsMeshParticleTracesField( *other.mParticleField ) )
  , mRendererContext( other.mRendererContext )
  , mFPS( other.mFPS )
  , mVpixMax( other.mVpixMax )
  , mParticleLifeTime( other.mParticleLifeTime )
 
{
 
}
 
 
void QgsMeshVectorTraceAnimationGenerator::seedRandomParticles( int count )
{
  mParticleField->setParticlesCount( count );
  mParticleField->addRandomParticles();
}
 
QImage QgsMeshVectorTraceAnimationGenerator::imageRendered()
{
  mParticleField->moveParticles();
  return mParticleField->image();
}
 
void QgsMeshVectorTraceAnimationGenerator::setFPS( int FPS )
{
  if ( FPS > 0 )
    mFPS = FPS;
  else
    mFPS = 1;
 
  updateFieldParameter();
}
 
void QgsMeshVectorTraceAnimationGenerator::setMaxSpeedPixel( int max )
{
  mVpixMax = max;
  updateFieldParameter();
}
 
void QgsMeshVectorTraceAnimationGenerator::setParticlesLifeTime( double particleLifeTime )
{
  mParticleLifeTime = particleLifeTime;
  updateFieldParameter();
}
 
void QgsMeshVectorTraceAnimationGenerator::setParticlesColor( const QColor &c )
{
  mParticleField->setParticlesColor( c );
}
 
void QgsMeshVectorTraceAnimationGenerator::setParticlesSize( double width )
{
  mParticleField->setParticleSize( width );
}
 
void QgsMeshVectorTraceAnimationGenerator::setTailFactor( double fct )
{
  mParticleField->setTailFactor( fct );
}
 
void QgsMeshVectorTraceAnimationGenerator::setMinimumTailLength( int l )
{
  mParticleField->setMinTailLength( l );
}
 
void QgsMeshVectorTraceAnimationGenerator::setTailPersitence( double p )
{
  if ( p < 0 )
    p = 0;
  if ( p > 1 )
    p = 1;
  mParticleField->setStumpFactor( int( 255 * p ) );
}
 
QgsMeshVectorTraceAnimationGenerator &QgsMeshVectorTraceAnimationGenerator::operator=( const QgsMeshVectorTraceAnimationGenerator &other )
{
  mParticleField.reset( new QgsMeshParticleTracesField( *( other.mParticleField ) ) );
  const_cast<QgsRenderContext &>( mRendererContext ) = other.mRendererContext;
  mFPS = other.mFPS;
  mVpixMax = other.mVpixMax;
  mParticleLifeTime = other.mParticleLifeTime;
 
  return ( *this );
}
 
void QgsMeshVectorTraceAnimationGenerator::updateFieldParameter()
{
  double fieldTimeStep = mVpixMax / static_cast<double>( mFPS );
  double fieldLifeTime = mParticleLifeTime * mFPS * fieldTimeStep;
  mParticleField->setTimeStep( fieldTimeStep );
  mParticleField->setParticlesLifeTime( fieldLifeTime );
}
 
QgsMeshVectorTraceRenderer::QgsMeshVectorTraceRenderer(
  const QgsTriangularMesh &triangularMesh,
  const QgsMeshDataBlock &dataSetVectorValues,
  const QgsMeshDataBlock &scalarActiveFaceFlagValues,
  bool dataIsOnVertices,
  const QgsMeshRendererVectorSettings &settings,
  QgsRenderContext &rendererContext,
  const QgsRectangle &layerExtent,
  double magMax )
  :  mParticleField( new QgsMeshParticleTracesField(
                       triangularMesh,
                       dataSetVectorValues,
                       scalarActiveFaceFlagValues,
                       layerExtent,
                       magMax,
                       dataIsOnVertices,
                       rendererContext,
                       settings.vectorStrokeColoring() ) )
  , mRendererContext( rendererContext )
{
 
  mParticleField->updateSize( rendererContext ) ;
 
  mParticleField->setParticleSize( rendererContext.convertToPainterUnits(
                                     settings.lineWidth(), Qgis::RenderUnit::Millimeters ) );
  mParticleField->setParticlesCount( settings.tracesSettings().particlesCount() );
  mParticleField->setTailFactor( 1 );
  mParticleField->setStumpParticleWithLifeTime( false );
  mParticleField->setTimeStep( rendererContext.convertToPainterUnits( settings.tracesSettings().maximumTailLength(),
                               settings.tracesSettings().maximumTailLengthUnit() ) ); //as the particles go through 1 pix for dt=1 and Vmax
  mParticleField->addRandomParticles();
  mParticleField->moveParticles();
}
 
void QgsMeshVectorTraceRenderer::draw()
{
  if ( mRendererContext.renderingStopped() )
    return;
  mRendererContext.painter()->drawImage( mParticleField->topLeft(), mParticleField->image() );
}