Branch data     Line data    Source code

       1                 :            : /***************************************************************************
       2                 :            :                               qgstininterpolator.cpp
       3                 :            :                               ----------------------
       4                 :            :   begin                : March 10, 2008
       5                 :            :   copyright            : (C) 2008 by Marco Hugentobler
       6                 :            :   email                : marco dot hugentobler at karto dot baug dot ethz dot ch
       7                 :            :  ***************************************************************************/
       8                 :            : 
       9                 :            : /***************************************************************************
      10                 :            :  *                                                                         *
      11                 :            :  *   This program is free software; you can redistribute it and/or modify  *
      12                 :            :  *   it under the terms of the GNU General Public License as published by  *
      13                 :            :  *   the Free Software Foundation; either version 2 of the License, or     *
      14                 :            :  *   (at your option) any later version.                                   *
      15                 :            :  *                                                                         *
      16                 :            :  ***************************************************************************/
      17                 :            : 
      18                 :            : #include "qgstininterpolator.h"
      19                 :            : #include "qgsfeatureiterator.h"
      20                 :            : #include "CloughTocherInterpolator.h"
      21                 :            : #include "qgsdualedgetriangulation.h"
      22                 :            : #include "NormVecDecorator.h"
      23                 :            : #include "LinTriangleInterpolator.h"
      24                 :            : #include "qgspoint.h"
      25                 :            : #include "qgsfeature.h"
      26                 :            : #include "qgsgeometry.h"
      27                 :            : #include "qgsvectorlayer.h"
      28                 :            : #include "qgswkbptr.h"
      29                 :            : #include "qgsfeedback.h"
      30                 :            : #include "qgscurve.h"
      31                 :            : #include "qgsmulticurve.h"
      32                 :            : #include "qgscurvepolygon.h"
      33                 :            : #include "qgsmultisurface.h"
      34                 :            : 
      35                 :          0 : QgsTinInterpolator::QgsTinInterpolator( const QList<LayerData> &inputData, TinInterpolation interpolation, QgsFeedback *feedback )
      36                 :          0 :   : QgsInterpolator( inputData )
      37                 :          0 :   , mIsInitialized( false )
      38                 :          0 :   , mFeedback( feedback )
      39                 :          0 :   , mInterpolation( interpolation )
      40                 :          0 : {
      41                 :          0 : }
      42                 :            : 
      43                 :          0 : QgsTinInterpolator::~QgsTinInterpolator()
      44                 :          0 : {
      45                 :          0 :   delete mTriangulation;
      46                 :          0 :   delete mTriangleInterpolator;
      47                 :          0 : }
      48                 :            : 
      49                 :          0 : int QgsTinInterpolator::interpolatePoint( double x, double y, double &result, QgsFeedback * )
      50                 :            : {
      51                 :          0 :   if ( !mIsInitialized )
      52                 :            :   {
      53                 :          0 :     initialize();
      54                 :          0 :   }
      55                 :            : 
      56                 :          0 :   if ( !mTriangleInterpolator )
      57                 :            :   {
      58                 :          0 :     return 1;
      59                 :            :   }
      60                 :            : 
      61                 :          0 :   QgsPoint r( 0, 0, 0 );
      62                 :          0 :   if ( !mTriangleInterpolator->calcPoint( x, y, r ) )
      63                 :            :   {
      64                 :          0 :     return 2;
      65                 :            :   }
      66                 :          0 :   result = r.z();
      67                 :          0 :   return 0;
      68                 :          0 : }
      69                 :            : 
      70                 :          0 : QgsFields QgsTinInterpolator::triangulationFields()
      71                 :            : {
      72                 :          0 :   return QgsTriangulation::triangulationFields();
      73                 :            : }
      74                 :            : 
      75                 :          0 : void QgsTinInterpolator::setTriangulationSink( QgsFeatureSink *sink )
      76                 :            : {
      77                 :          0 :   mTriangulationSink = sink;
      78                 :          0 : }
      79                 :          0 : 
      80                 :          0 : void QgsTinInterpolator::initialize()
      81                 :            : {
      82                 :          0 :   QgsDualEdgeTriangulation *dualEdgeTriangulation = new QgsDualEdgeTriangulation( 100000 );
      83                 :          0 :   if ( mInterpolation == CloughTocher )
      84                 :            :   {
      85                 :          0 :     NormVecDecorator *dec = new NormVecDecorator();
      86                 :          0 :     dec->addTriangulation( dualEdgeTriangulation );
      87                 :          0 :     mTriangulation = dec;
      88                 :          0 :   }
      89                 :            :   else
      90                 :            :   {
      91                 :          0 :     mTriangulation = dualEdgeTriangulation;
      92                 :            :   }
      93                 :            : 
      94                 :            :   //get number of features if we use a progress bar
      95                 :          0 :   int nFeatures = 0;
      96                 :          0 :   int nProcessedFeatures = 0;
      97                 :          0 :   if ( mFeedback )
      98                 :            :   {
      99                 :          0 :     for ( const LayerData &layer :  std::as_const( mLayerData ) )
     100                 :            :     {
     101                 :          0 :       if ( layer.source )
     102                 :            :       {
     103                 :          0 :         nFeatures += layer.source->featureCount();
     104                 :          0 :       }
     105                 :            :     }
     106                 :          0 :   }
     107                 :            : 
     108                 :          0 :   const QgsCoordinateReferenceSystem crs = !mLayerData.empty() ? mLayerData.at( 0 ).source->sourceCrs() : QgsCoordinateReferenceSystem();
     109                 :            : 
     110                 :          0 :   QgsFeature f;
     111                 :          0 :   for ( const LayerData &layer : std::as_const( mLayerData ) )
     112                 :            :   {
     113                 :          0 :     if ( layer.source )
     114                 :            :     {
     115                 :          0 :       QgsAttributeList attList;
     116                 :          0 :       switch ( layer.valueSource )
     117                 :            :       {
     118                 :            :         case QgsInterpolator::ValueAttribute:
     119                 :          0 :           attList.push_back( layer.interpolationAttribute );
     120                 :          0 :           break;
     121                 :            : 
     122                 :            :         case QgsInterpolator::ValueM:
     123                 :            :         case QgsInterpolator::ValueZ:
     124                 :          0 :           break;
     125                 :            :       }
     126                 :            : 
     127                 :          0 :       QgsFeatureIterator fit = layer.source->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( attList ).setDestinationCrs( crs, layer.transformContext ) );
     128                 :            : 
     129                 :          0 :       while ( fit.nextFeature( f ) )
     130                 :            :       {
     131                 :          0 :         if ( mFeedback )
     132                 :            :         {
     133                 :          0 :           if ( mFeedback->isCanceled() )
     134                 :            :           {
     135                 :          0 :             break;
     136                 :            :           }
     137                 :          0 :           if ( nFeatures > 0 )
     138                 :          0 :             mFeedback->setProgress( 100.0 * static_cast< double >( nProcessedFeatures ) / nFeatures );
     139                 :          0 :         }
     140                 :          0 :         insertData( f, layer.valueSource, layer.interpolationAttribute, layer.sourceType );
     141                 :          0 :         ++nProcessedFeatures;
     142                 :            :       }
     143                 :          0 :     }
     144                 :            :   }
     145                 :            : 
     146                 :          0 :   if ( mInterpolation == CloughTocher )
     147                 :            :   {
     148                 :          0 :     CloughTocherInterpolator *ctInterpolator = new CloughTocherInterpolator();
     149                 :          0 :     NormVecDecorator *dec = dynamic_cast<NormVecDecorator *>( mTriangulation );
     150                 :          0 :     if ( dec )
     151                 :            :     {
     152                 :          0 :       dec->estimateFirstDerivatives( mFeedback );
     153                 :          0 :       ctInterpolator->setTriangulation( dec );
     154                 :          0 :       dec->setTriangleInterpolator( ctInterpolator );
     155                 :          0 :       mTriangleInterpolator = ctInterpolator;
     156                 :          0 :     }
     157                 :          0 :   }
     158                 :            :   else //linear
     159                 :            :   {
     160                 :          0 :     mTriangleInterpolator = new LinTriangleInterpolator( dualEdgeTriangulation );
     161                 :            :   }
     162                 :          0 :   mIsInitialized = true;
     163                 :            : 
     164                 :            :   //debug
     165                 :          0 :   if ( mTriangulationSink )
     166                 :            :   {
     167                 :          0 :     dualEdgeTriangulation->saveTriangulation( mTriangulationSink, mFeedback );
     168                 :          0 :   }
     169                 :          0 : }
     170                 :            : 
     171                 :          0 : int QgsTinInterpolator::insertData( const QgsFeature &f, QgsInterpolator::ValueSource source, int attr, SourceType type )
     172                 :            : {
     173                 :          0 :   QgsGeometry g = f.geometry();
     174                 :          0 :   if ( g.isNull() || g.isEmpty() )
     175                 :            :   {
     176                 :          0 :     return 2;
     177                 :            :   }
     178                 :            : 
     179                 :            :   //check attribute value
     180                 :          0 :   double attributeValue = 0;
     181                 :          0 :   bool attributeConversionOk = false;
     182                 :          0 :   switch ( source )
     183                 :            :   {
     184                 :            :     case ValueAttribute:
     185                 :            :     {
     186                 :          0 :       QVariant attributeVariant = f.attribute( attr );
     187                 :          0 :       if ( !attributeVariant.isValid() ) //attribute not found, something must be wrong (e.g. NULL value)
     188                 :            :       {
     189                 :          0 :         return 3;
     190                 :            :       }
     191                 :          0 :       attributeValue = attributeVariant.toDouble( &attributeConversionOk );
     192                 :          0 :       if ( !attributeConversionOk || std::isnan( attributeValue ) ) //don't consider vertices with attributes like 'nan' for the interpolation
     193                 :            :       {
     194                 :          0 :         return 4;
     195                 :            :       }
     196                 :          0 :       break;
     197                 :          0 :     }
     198                 :            : 
     199                 :            :     case ValueM:
     200                 :          0 :       if ( !g.constGet()->isMeasure() )
     201                 :          0 :         return 3;
     202                 :            :       else
     203                 :          0 :         break;
     204                 :            : 
     205                 :            :     case ValueZ:
     206                 :          0 :       if ( !g.constGet()->is3D() )
     207                 :          0 :         return 3;
     208                 :            :       else
     209                 :          0 :         break;
     210                 :            :   }
     211                 :            : 
     212                 :            : 
     213                 :          0 :   switch ( type )
     214                 :            :   {
     215                 :            :     case SourcePoints:
     216                 :            :     {
     217                 :          0 :       if ( addPointsFromGeometry( g, source, attributeValue ) != 0 )
     218                 :          0 :         return -1;
     219                 :          0 :       break;
     220                 :            :     }
     221                 :            : 
     222                 :            :     case SourceBreakLines:
     223                 :            :     case SourceStructureLines:
     224                 :            :     {
     225                 :          0 :       switch ( QgsWkbTypes::geometryType( g.wkbType() ) )
     226                 :            :       {
     227                 :            :         case QgsWkbTypes::PointGeometry:
     228                 :            :         {
     229                 :          0 :           if ( addPointsFromGeometry( g, source, attributeValue ) != 0 )
     230                 :          0 :             return -1;
     231                 :          0 :           break;
     232                 :            :         }
     233                 :            : 
     234                 :            :         case QgsWkbTypes::LineGeometry:
     235                 :            :         case QgsWkbTypes::PolygonGeometry:
     236                 :            :         {
     237                 :            :           // need to extract all rings from input geometry
     238                 :          0 :           std::vector<const QgsCurve *> curves;
     239                 :          0 :           if ( QgsWkbTypes::geometryType( g.wkbType() ) == QgsWkbTypes::PolygonGeometry )
     240                 :            :           {
     241                 :          0 :             std::vector< const QgsCurvePolygon * > polygons;
     242                 :          0 :             if ( g.isMultipart() )
     243                 :            :             {
     244                 :          0 :               const QgsMultiSurface *ms = qgsgeometry_cast< const QgsMultiSurface * >( g.constGet() );
     245                 :          0 :               for ( int i = 0; i < ms->numGeometries(); ++i )
     246                 :            :               {
     247                 :          0 :                 polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( ms->geometryN( i ) ) );
     248                 :          0 :               }
     249                 :          0 :             }
     250                 :            :             else
     251                 :            :             {
     252                 :          0 :               polygons.emplace_back( qgsgeometry_cast< const QgsCurvePolygon * >( g.constGet() ) );
     253                 :            :             }
     254                 :            : 
     255                 :          0 :             for ( const QgsCurvePolygon *polygon : polygons )
     256                 :            :             {
     257                 :          0 :               if ( !polygon )
     258                 :          0 :                 continue;
     259                 :            : 
     260                 :          0 :               if ( polygon->exteriorRing() )
     261                 :          0 :                 curves.emplace_back( polygon->exteriorRing() );
     262                 :            : 
     263                 :          0 :               for ( int i = 0; i < polygon->numInteriorRings(); ++i )
     264                 :            :               {
     265                 :          0 :                 curves.emplace_back( polygon->interiorRing( i ) );
     266                 :          0 :               }
     267                 :            :             }
     268                 :          0 :           }
     269                 :            :           else
     270                 :            :           {
     271                 :          0 :             if ( g.isMultipart() )
     272                 :            :             {
     273                 :          0 :               const QgsMultiCurve *mc = qgsgeometry_cast< const QgsMultiCurve * >( g.constGet() );
     274                 :          0 :               for ( int i = 0; i < mc->numGeometries(); ++i )
     275                 :            :               {
     276                 :          0 :                 curves.emplace_back( mc->curveN( i ) );
     277                 :          0 :               }
     278                 :          0 :             }
     279                 :            :             else
     280                 :            :             {
     281                 :          0 :               curves.emplace_back( qgsgeometry_cast< const QgsCurve * >( g.constGet() ) );
     282                 :            :             }
     283                 :            :           }
     284                 :            : 
     285                 :          0 :           for ( const QgsCurve *curve : curves )
     286                 :            :           {
     287                 :          0 :             if ( !curve )
     288                 :          0 :               continue;
     289                 :            : 
     290                 :          0 :             QgsPointSequence linePoints;
     291                 :          0 :             curve->points( linePoints );
     292                 :          0 :             for ( QgsPoint &point : linePoints )
     293                 :            :             {
     294                 :          0 :               switch ( source )
     295                 :            :               {
     296                 :            :                 case ValueAttribute:
     297                 :          0 :                   if ( point.is3D() )
     298                 :          0 :                     point.setZ( attributeValue );
     299                 :            :                   else
     300                 :          0 :                     point.addZValue( attributeValue );
     301                 :          0 :                   break;
     302                 :            : 
     303                 :            :                 case ValueM:
     304                 :          0 :                   if ( point.is3D() )
     305                 :          0 :                     point.setZ( point.m() );
     306                 :            :                   else
     307                 :          0 :                     point.addZValue( point.m() );
     308                 :          0 :                   break;
     309                 :            : 
     310                 :            :                 case ValueZ:
     311                 :          0 :                   break;
     312                 :            :               }
     313                 :            :             }
     314                 :          0 :             mTriangulation->addLine( linePoints, type );
     315                 :          0 :           }
     316                 :            :           break;
     317                 :          0 :         }
     318                 :            :         case QgsWkbTypes::UnknownGeometry:
     319                 :            :         case QgsWkbTypes::NullGeometry:
     320                 :          0 :           break;
     321                 :            :       }
     322                 :          0 :       break;
     323                 :            :     }
     324                 :            :   }
     325                 :            : 
     326                 :          0 :   return 0;
     327                 :          0 : }
     328                 :            : 
     329                 :            : 
     330                 :          0 : int QgsTinInterpolator::addPointsFromGeometry( const QgsGeometry &g, ValueSource source, double attributeValue )
     331                 :            : {
     332                 :            :   // loop through all vertices and add to triangulation
     333                 :          0 :   for ( auto point = g.vertices_begin(); point != g.vertices_end(); ++point )
     334                 :            :   {
     335                 :          0 :     QgsPoint p = *point;
     336                 :          0 :     double z = 0;
     337                 :          0 :     switch ( source )
     338                 :            :     {
     339                 :            :       case ValueAttribute:
     340                 :          0 :         z = attributeValue;
     341                 :          0 :         break;
     342                 :            : 
     343                 :            :       case ValueZ:
     344                 :          0 :         z = p.z();
     345                 :          0 :         break;
     346                 :            : 
     347                 :            :       case ValueM:
     348                 :          0 :         z = p.m();
     349                 :          0 :         break;
     350                 :            :     }
     351                 :          0 :     if ( mTriangulation->addPoint( QgsPoint( p.x(), p.y(), z ) ) == -100 )
     352                 :            :     {
     353                 :          0 :       return -1;
     354                 :            :     }
     355                 :          0 :   }
     356                 :          0 :   return 0;
     357                 :          0 : }