Fixed few bugs + minor rendering improvements

- using a color ramp to render points - only rendering nodes that intersect the map extent
qgis · Oct 26, 2020 · da969e8 · da969e8
1 parent 1d4b8e7
commit da969e8
Show file tree

Hide file tree

Showing 7 changed files with 108 additions and 45 deletions.
diff --git a/python/core/auto_generated/pointcloud/qgspointcloudrenderer.sip.in b/python/core/auto_generated/pointcloud/qgspointcloudrenderer.sip.in
@@ -12,6 +12,8 @@
 
 
 
+
+
 class QgsPointCloudRenderer: QgsMapLayerRenderer
 {
 %Docstring

diff --git a/src/core/pointcloud/qgspointclouddecoder.cpp b/src/core/pointcloud/qgspointclouddecoder.cpp
@@ -138,4 +138,5 @@ QVector<qint32> QgsPointCloudDecoder::decompressLaz( const QString &filename )
   }
   float t = common::since( start );
   std::cout << "LAZ-PERF Read through the points in " << t << " seconds." << std::endl;
+  return data;
 }
diff --git a/src/core/pointcloud/qgspointcloudindex.cpp b/src/core/pointcloud/qgspointcloudindex.cpp
@@ -53,12 +53,12 @@ uint qHash( const IndexedPointCloudNode &id )
 
 QgsPointCloudDataBounds::QgsPointCloudDataBounds() = default;
 
-QgsPointCloudDataBounds::QgsPointCloudDataBounds( qint32 xmin, qint32 ymin, qint32 xmax, qint32 ymax, qint32 zmin, qint32 zmax )
+QgsPointCloudDataBounds::QgsPointCloudDataBounds( qint32 xmin, qint32 ymin, qint32 zmin, qint32 xmax, qint32 ymax, qint32 zmax )
   : mXMin( xmin )
   , mYMin( ymin )
+  , mZMin( zmin )
   , mXMax( xmax )
   , mYMax( ymax )
-  , mZMin( zmin )
   , mZMax( zmax )
 {
 
@@ -67,9 +67,9 @@ QgsPointCloudDataBounds::QgsPointCloudDataBounds( qint32 xmin, qint32 ymin, qint
 QgsPointCloudDataBounds::QgsPointCloudDataBounds( const QgsPointCloudDataBounds &obj )
   : mXMin( obj.xMin() )
   , mYMin( obj.yMin() )
+  , mZMin( obj.zMin() )
   , mXMax( obj.xMax() )
   , mYMax( obj.yMax() )
-  , mZMin( obj.zMin() )
   , mZMax( obj.zMax() )
 {
 }
@@ -104,6 +104,14 @@ qint32 QgsPointCloudDataBounds::zMax() const
   return mZMax;
 }
 
+QgsRectangle QgsPointCloudDataBounds::mapExtent( const QgsVector3D &offset, const QgsVector3D &scale ) const
+{
+  return QgsRectangle(
+           mXMin * scale.x() + offset.x(), mYMin * scale.y() + offset.y(),
+           mXMax * scale.x() + offset.x(), mYMax * scale.y() + offset.y()
+         );
+}
+
 QgsPointCloudIndex::QgsPointCloudIndex() = default;
 
 QgsPointCloudIndex::~QgsPointCloudIndex() = default;
@@ -146,6 +154,10 @@ bool QgsPointCloudIndex::load( const QString &fileName )
   QJsonArray bounds_conforming = doc["boundsConforming"].toArray();
   if ( bounds.size() != 6 )
     return false;
+  mExtent.set( bounds_conforming[0].toDouble(), bounds_conforming[1].toDouble(),
+               bounds_conforming[3].toDouble(), bounds_conforming[4].toDouble() );
+  mZMin = bounds_conforming[2].toDouble();
+  mZMax = bounds_conforming[5].toDouble();
 
   QJsonArray schemaArray = doc["schema"].toArray();
 
@@ -195,10 +207,10 @@ bool QgsPointCloudIndex::load( const QString &fileName )
 
   mRootBounds = QgsPointCloudDataBounds(
                   ( xmin - mOffset.x() ) / mScale.x(),
-                  ( xmax - mOffset.x() ) / mScale.x(),
                   ( ymin - mOffset.y() ) / mScale.y(),
-                  ( ymax - mOffset.y() ) / mScale.y(),
                   ( zmin - mOffset.z() ) / mScale.z(),
+                  ( xmax - mOffset.x() ) / mScale.x(),
+                  ( ymax - mOffset.y() ) / mScale.y(),
                   ( zmax - mOffset.z() ) / mScale.z()
                 );
 
@@ -267,13 +279,17 @@ QVector<qint32> QgsPointCloudIndex::nodePositionDataAsInt32( const IndexedPointC
   {
     QString filename = QString( "%1/ept-data/%2.zst" ).arg( mDirectory ).arg( n.toString() );
     Q_ASSERT( QFile::exists( filename ) );
-    return QgsPointCloudDecoder::decompressBinary( filename );
+    return QgsPointCloudDecoder::decompressZStandard( filename );
   }
-  else     // if ( mDataType == "laz" )
+  else if ( mDataType == "laszip" )
   {
     QString filename = QString( "%1/ept-data/%2.laz" ).arg( mDirectory ).arg( n.toString() );
     Q_ASSERT( QFile::exists( filename ) );
-    return QgsPointCloudDecoder::decompressBinary( filename );
+    return QgsPointCloudDecoder::decompressLaz( filename );
+  }
+  else
+  {
+    Q_ASSERT( false );  // unsupported
   }
 }
 
@@ -292,10 +308,15 @@ QgsPointCloudDataBounds QgsPointCloudIndex::nodeBounds( const IndexedPointCloudN
   zMin = round( mRootBounds.zMin() + dLevel * n.z );
   zMax = round( mRootBounds.zMin() + dLevel * ( n.z + 1 ) );
 
-  QgsPointCloudDataBounds db( xMin, xMax, yMin, yMax, zMin, zMax );
+  QgsPointCloudDataBounds db( xMin, yMin, zMin, xMax, yMax, zMax );
   return db;
 }
 
+QgsRectangle QgsPointCloudIndex::nodeMapExtent( const IndexedPointCloudNode &n )
+{
+  return nodeBounds( n ).mapExtent( mOffset, mScale );
+}
+
 QString QgsPointCloudIndex::wkt() const
 {
   return mWkt;

diff --git a/src/core/pointcloud/qgspointcloudindex.h b/src/core/pointcloud/qgspointcloudindex.h
@@ -26,6 +26,7 @@
 #include <QList>
 
 #include "qgis_core.h"
+#include "qgsrectangle.h"
 #include "qgsvector3d.h"
 #include "qgis_sip.h"
 
@@ -66,23 +67,26 @@ class CORE_EXPORT QgsPointCloudDataBounds
 {
   public:
     QgsPointCloudDataBounds(); // invalid
-    QgsPointCloudDataBounds( qint32 xmin, qint32 ymin, qint32 xmax, qint32 ymax, qint32 zmin, qint32 zmax );
+    QgsPointCloudDataBounds( qint32 xmin, qint32 ymin, qint32 zmin, qint32 xmax, qint32 ymax, qint32 zmax );
     QgsPointCloudDataBounds( const QgsPointCloudDataBounds &obj );
 
     qint32 xMin() const;
 
     qint32 yMin() const;
 
+    qint32 zMin() const;
+
     qint32 xMax() const;
 
     qint32 yMax() const;
 
-    qint32 zMin() const;
-
     qint32 zMax() const;
 
+    //! Returns 2D rectangle in map coordinates
+    QgsRectangle mapExtent( const QgsVector3D &offset, const QgsVector3D &scale ) const;
+
   private:
-    qint32 mXMin, mYMin, mXMax, mYMax, mZMin, mZMax;
+    qint32 mXMin, mYMin, mZMin, mXMax, mYMax, mZMax;
 };
 
 /**
@@ -109,7 +113,11 @@ class CORE_EXPORT QgsPointCloudIndex: public QObject
 
     QVector<qint32> nodePositionDataAsInt32( const IndexedPointCloudNode &n );
 
+    QgsRectangle extent() const { return mExtent; }
+    double zMin() const { return mZMin; }
+    double zMax() const { return mZMax; }
     QgsPointCloudDataBounds nodeBounds( const IndexedPointCloudNode &n );
+    QgsRectangle nodeMapExtent( const IndexedPointCloudNode &n );
     QString wkt() const;
 
     QgsVector3D scale() const;
@@ -120,6 +128,8 @@ class CORE_EXPORT QgsPointCloudIndex: public QObject
     QString mDirectory;
     QString mDataType;
 
+    QgsRectangle mExtent;  //!< 2D extent of data
+    double mZMin = 0, mZMax = 0;   //!< Vertical extent of data
     QHash<IndexedPointCloudNode, int> mHierarchy;
     QgsVector3D mScale; //!< Scale of our int32 coordinates compared to CRS coords
     QgsVector3D mOffset; //!< Offset of our int32 coordinates compared to CRS coords

diff --git a/src/core/pointcloud/qgspointcloudlayer.cpp b/src/core/pointcloud/qgspointcloudlayer.cpp
@@ -40,16 +40,7 @@ QgsRectangle QgsPointCloudLayer::extent() const
   if ( !mPointCloudIndex )
     return QgsRectangle();
 
-  QgsPointCloudDataBounds bounds = mPointCloudIndex->nodeBounds( mPointCloudIndex->root() );
-  QgsVector3D offset = mPointCloudIndex->offset();
-  QgsVector3D scale = mPointCloudIndex->scale();
-
-  double xMin = offset.x() + scale.x() * bounds.xMin();
-  double yMin = offset.y() + scale.y() * bounds.yMin();
-  double xMax = offset.x() + scale.x() * bounds.xMax();
-  double yMax = offset.y() + scale.y() * bounds.yMax();
-
-  return QgsRectangle( xMin, yMin, xMax, yMax );
+  return mPointCloudIndex->extent();
 }
 
 QgsMapLayerRenderer *QgsPointCloudLayer::createMapRenderer( QgsRenderContext &rendererContext )

diff --git a/src/core/pointcloud/qgspointcloudrenderer.cpp b/src/core/pointcloud/qgspointcloudrenderer.cpp
@@ -21,23 +21,36 @@
 #include "qgspointcloudlayer.h"
 #include "qgsrendercontext.h"
 #include "qgspointcloudindex.h"
+#include "qgsstyle.h"
+#include "qgscolorramp.h"
 
 QgsPointCloudRenderer::QgsPointCloudRenderer( QgsPointCloudLayer *layer, QgsRenderContext &context )
   : QgsMapLayerRenderer( layer->id(), &context )
   , mLayer( layer )
 {
-
+  // TODO: we must not keep pointer to mLayer (it's dangerous) - we must copy anything we need for rendering
+  // or use some locking to prevent read/write from multiple threads
+
+  // TODO: use config from layer
+  mConfig.penWidth = context.convertToPainterUnits( 1, QgsUnitTypes::RenderUnit::RenderMillimeters );
+  mConfig.zMin = 400;
+  mConfig.zMax = 600;
+  mConfig.colorRamp.reset( QgsStyle::defaultStyle()->colorRamp( "Viridis" ) );
 }
 
-static QList<IndexedPointCloudNode> _traverseTree( QgsPointCloudIndex *pc, IndexedPointCloudNode n, int maxDepth )
+static QList<IndexedPointCloudNode> _traverseTree( QgsPointCloudIndex *pc, const QgsRectangle &extent, IndexedPointCloudNode n, int maxDepth )
 {
   QList<IndexedPointCloudNode> nodes;
+
+  if ( !extent.intersects( pc->nodeMapExtent( n ) ) )
+    return nodes;
+
   nodes.append( n );
 
   for ( auto nn : pc->children( n ) )
   {
     if ( maxDepth > 0 )
-      nodes += _traverseTree( pc, nn, maxDepth - 1 );
+      nodes += _traverseTree( pc, extent, nn, maxDepth - 1 );
   }
 
   return nodes;
@@ -56,27 +69,24 @@ bool QgsPointCloudRenderer::render()
   painter->save();
   context.setPainterFlagsUsingContext( painter );
 
-  QPen pen = painter->pen();
-  pen.setCapStyle( Qt::FlatCap );
-  pen.setJoinStyle( Qt::MiterJoin );
-
-  double penWidth = context.convertToPainterUnits( 1, QgsUnitTypes::RenderUnit::RenderMillimeters );
-  pen.setWidthF( penWidth );
-  pen.setColor( Qt::black );
-  painter->setPen( pen );
-
   QgsPointCloudDataBounds db;
 
   QElapsedTimer t;
   t.start();
 
-  // TODO: traverse with spatial filter
-  QList<IndexedPointCloudNode> nodes = _traverseTree( pc, pc->root(), 4 );
+  // TODO: set depth based on map units per pixel
+  int depth = 3;
+  QList<IndexedPointCloudNode> nodes = _traverseTree( pc, context.mapExtent(), pc->root(), depth );
 
   // drawing
   for ( auto n : nodes )
   {
-    drawData( painter, pc->nodePositionDataAsInt32( n ) );
+    if ( context.renderingStopped() )
+    {
+      qDebug() << "canceled";
+      break;
+    }
+    drawData( painter, pc->nodePositionDataAsInt32( n ), mConfig );
   }
 
   qDebug() << "totals:" << nodesDrawn << "nodes | " << pointsDrawn << " points | " << t.elapsed() << "ms";
@@ -100,26 +110,38 @@ void QgsPointCloudRenderer::readXml( const QDomElement &elem, const QgsReadWrite
 
 }
 
-void QgsPointCloudRenderer::drawData( QPainter *painter, const QVector<qint32> &data )
+void QgsPointCloudRenderer::drawData( QPainter *painter, const QVector<qint32> &data, const QgsPointCloudRendererConfig &config )
 {
   const QgsMapToPixel mapToPixel = renderContext()->mapToPixel();
   const QgsVector3D scale = mLayer->pointCloudIndex()->scale();
   const QgsVector3D offset = mLayer->pointCloudIndex()->offset();
 
+  QgsRectangle mapExtent = renderContext()->mapExtent();
+
+  QPen pen;
+  pen.setWidth( config.penWidth );
+  pen.setCapStyle( Qt::FlatCap );
+  //pen.setJoinStyle( Qt::MiterJoin );
+
   const qint32 *ptr = data.constData();
   int count = data.count() / 3;
   for ( int i = 0; i < count; ++i )
   {
     qint32 ix = ptr[i * 3 + 0];
     qint32 iy = ptr[i * 3 + 1];
-    // qint32 iz = ptr[i * 3 + 2];
+    qint32 iz = ptr[i * 3 + 2];
 
     double x = offset.x() + scale.x() * ix;
     double y = offset.y() + scale.y() * iy;
-    // double z = offset.z() + scale.z() * iz;
-    mapToPixel.transformInPlace( x, y );
-
-    painter->drawPoint( QPointF( x, y ) );
+    if ( mapExtent.contains( QgsPointXY( x, y ) ) )
+    {
+      double z = offset.z() + scale.z() * iz;
+      mapToPixel.transformInPlace( x, y );
+
+      pen.setColor( config.colorRamp->color( ( z - config.zMin ) / ( config.zMax - config.zMin ) ) );
+      painter->setPen( pen );
+      painter->drawPoint( QPointF( x, y ) );
+    }
   }
 
   // stats

diff --git a/src/core/pointcloud/qgspointcloudrenderer.h b/src/core/pointcloud/qgspointcloudrenderer.h
@@ -19,6 +19,7 @@
 #define QGSPOINTCLOUDRENDERER_H
 
 #include "qgis_core.h"
+#include "qgscolorramp.h"
 #include "qgsmaplayerrenderer.h"
 #include "qgsreadwritecontext.h"
 #include "qgspointcloudindex.h"
@@ -27,9 +28,22 @@
 #include <QString>
 #include <QPainter>
 
+
 class QgsRenderContext;
 class QgsPointCloudLayer;
 
+#ifndef SIP_RUN
+
+class QgsPointCloudRendererConfig
+{
+  public:
+    double zMin = 0, zMax = 0;
+    int penWidth = 1;
+    std::unique_ptr<QgsColorRamp> colorRamp;
+};
+
+#endif
+
 
 /**
  * \ingroup core
@@ -55,14 +69,16 @@ class CORE_EXPORT QgsPointCloudRenderer: public QgsMapLayerRenderer
   private:
     QgsPointCloudLayer *mLayer = nullptr;
 
+    QgsPointCloudRendererConfig mConfig;
+
     // int imgW, imgH; // DO WE NEED AT ALL?
     // QgsPointCloudDataBounds mBounds; // DO WE NEED AT ALL?
 
     // some stats
     int nodesDrawn = 0;
     int pointsDrawn = 0;
 
-    void drawData( QPainter *painter, const QVector<qint32> &data );
+    void drawData( QPainter *painter, const QVector<qint32> &data, const QgsPointCloudRendererConfig &config );
 };