port the basic renderer from Martin's prototype

qgis · Oct 26, 2020 · 8a42c57 · 8a42c57
1 parent 13ecb8c
commit 8a42c57
Show file tree

Hide file tree

Showing 16 changed files with 766 additions and 53 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -350,6 +350,9 @@ IF(WITH_CORE)
   FIND_PACKAGE(ZLIB REQUIRED)       # for decompression of vector tiles in MBTiles file
   MESSAGE(STATUS "Found zlib: ${ZLIB_LIBRARIES}")
 
+  FIND_PACKAGE(ZSTD REQUIRED) # for decompression of point clouds
+  FIND_PACKAGE(LazPerf REQUIRED) # for decompression of point clouds
+
   # optional
   IF (WITH_POSTGRESQL)
     FIND_PACKAGE(Postgres) # PostgreSQL provider

diff --git a/cmake/FindLazPerf.cmake b/cmake/FindLazPerf.cmake
@@ -0,0 +1,30 @@
+# CMake module to search for laz-perf
+#
+# Once done this will define
+#
+#  LazPerf_FOUND - system has the zip library
+#  LazPerf_INCLUDE_DIRS - the zip include directories
+#
+# Copyright (c) 2020, Peter Petrik, <zilolv@gmail.com>
+#
+# Redistribution and use is allowed according to the terms of the BSD license.
+# For details see the accompanying COPYING-CMAKE-SCRIPTS file.
+
+FIND_PATH(LazPerf_INCLUDE_DIR
+  laz-perf/io.hpp
+  "$ENV{LIB_DIR}/include"
+  "$ENV{INCLUDE}"
+  /usr/local/include
+  /usr/include
+)
+
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(LazPerf DEFAULT_MSG LazPerf_INCLUDE_DIR)
+
+MARK_AS_ADVANCED(LazPerf_INCLUDE_DIR)
+
+IF (LazPerf_FOUND)
+  MESSAGE(STATUS "Found laz-perf: ${LazPerf_INCLUDE_DIR}")
+ELSE (LazPerf_FOUND)
+  MESSAGE(FATAL_ERROR "Could not find laz-perf")
+ENDIF (LazPerf_FOUND)
diff --git a/cmake/FindZSTD.cmake b/cmake/FindZSTD.cmake
@@ -0,0 +1,34 @@
+# CMake module to search for libzstd
+#
+# Once done this will define
+#
+#  ZSTD_FOUND - system has the zip library
+#  ZSTD_INCLUDE_DIRS - the zip include directories
+#  ZSTD_LIBRARY - Link this to use the zip library
+#
+# Copyright (c) 2020, Peter Petrik, <zilolv@gmail.com>
+#
+# Redistribution and use is allowed according to the terms of the BSD license.
+# For details see the accompanying COPYING-CMAKE-SCRIPTS file.
+
+FIND_PATH(ZSTD_INCLUDE_DIR
+  zstd.h
+  "$ENV{LIB_DIR}/include"
+  "$ENV{INCLUDE}"
+  /usr/local/include
+  /usr/include
+)
+
+FIND_LIBRARY(ZSTD_LIBRARY NAMES zstd PATHS "$ENV{LIB_DIR}/lib" "$ENV{LIB}" /usr/local/lib /usr/lib )
+
+INCLUDE(FindPackageHandleStandardArgs)
+FIND_PACKAGE_HANDLE_STANDARD_ARGS(ZSTD DEFAULT_MSG
+                                  ZSTD_LIBRARY ZSTD_INCLUDE_DIR)
+
+MARK_AS_ADVANCED(ZSTD_LIBRARY ZSTD_INCLUDE_DIR)
+
+IF (ZSTD_FOUND)
+  MESSAGE(STATUS "Found ZSTD: ${ZSTD_LIBRARY}")
+ELSE (ZSTD_FOUND)
+  MESSAGE(FATAL_ERROR "Could not find ZSTD")
+ENDIF (ZSTD_FOUND)
diff --git a/python/core/auto_generated/pointcloud/qgspointcloudindex.sip.in b/python/core/auto_generated/pointcloud/qgspointcloudindex.sip.in
diff --git a/python/core/auto_generated/pointcloud/qgspointcloudrenderer.sip.in b/python/core/auto_generated/pointcloud/qgspointcloudrenderer.sip.in
@@ -11,6 +11,7 @@
 
 
 
+
 class QgsPointCloudRenderer: QgsMapLayerRenderer
 {
 %Docstring
@@ -38,7 +39,6 @@ Represents a 2D renderer of point cloud data
     void writeXml( QDomElement &elem, const QgsReadWriteContext &context ) const;
     void readXml( const QDomElement &elem, const QgsReadWriteContext &context );
 
-  protected:
 };
 
 

diff --git a/python/core/core_auto.sip b/python/core/core_auto.sip
@@ -439,7 +439,6 @@
 %Include auto_generated/mesh/qgsmeshcalculator.sip
 %Include auto_generated/pointcloud/qgspointcloudlayer.sip
 %Include auto_generated/pointcloud/qgspointcloudrenderer.sip
-%Include auto_generated/pointcloud/qgspointcloudindex.sip
 %Include auto_generated/metadata/qgsabstractmetadatabase.sip
 %Include auto_generated/metadata/qgslayermetadata.sip
 %Include auto_generated/metadata/qgslayermetadataformatter.sip

diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
@@ -626,6 +626,7 @@ SET(QGIS_CORE_SRCS
   pointcloud/qgspointcloudindex.cpp
   pointcloud/qgspointclouddataitems.cpp
   pointcloud/qgspointcloudprovidermetadata.cpp
+  pointcloud/qgspointclouddecoder.cpp
 
   labeling/qgslabelfeature.cpp
   labeling/qgslabelingengine.cpp
@@ -1301,6 +1302,7 @@ SET(QGIS_CORE_HDRS
   pointcloud/qgspointcloudindex.h
   pointcloud/qgspointclouddataitems.h
   pointcloud/qgspointcloudprovidermetadata.h
+  pointcloud/qgspointclouddecoder.h
 
   metadata/qgsabstractmetadatabase.h
   metadata/qgslayermetadata.h
@@ -1612,6 +1614,8 @@ INCLUDE_DIRECTORIES(SYSTEM
   ${Qt5SerialPort_INCLUDE_DIRS}
   ${Protobuf_INCLUDE_DIRS}
   ${ZLIB_INCLUDE_DIRS}
+  ${ZSTD_INCLUDE_DIR}
+  ${LazPerf_INCLUDE_DIR}
 )
 
 
@@ -1748,6 +1752,7 @@ TARGET_LINK_LIBRARIES(qgis_core
   ${SQLITE3_LIBRARY}
   ${SPATIALITE_LIBRARY}
   ${LIBZIP_LIBRARY}
+  ${ZSTD_LIBRARY}
   ${Protobuf_LITE_LIBRARY}
   ${ZLIB_LIBRARIES}
 )

diff --git a/src/core/pointcloud/qgspointclouddecoder.cpp b/src/core/pointcloud/qgspointclouddecoder.cpp
@@ -0,0 +1,169 @@
+/***************************************************************************
+                         qgspointcloudrenderer.cpp
+                         --------------------
+    begin                : October 2020
+    copyright            : (C) 2020 by Peter Petrik
+    email                : zilolv 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 "qgspointclouddecoder.h"
+#include "qgspointcloudindex.h"
+#include "qgsvector3d.h"
+
+#include <zstd.h>
+#include <QFile>
+#include <iostream>
+
+#include "laz-perf/io.hpp"
+#include "laz-perf/common/common.hpp"
+
+QVector<qint32> QgsPointCloudDecoder::decompressBinary(const QString& filename, QgsPointCloudDataBounds &db)
+{
+  Q_ASSERT( QFile::exists( filename ) );
+
+  QFile f( filename );
+  bool r = f.open(QIODevice::ReadOnly);
+  Q_ASSERT(r);
+
+  // WHY??? per-record should be 18 based on schema, not 46
+  int stride = 46; //18;
+  int count = f.size() / stride;
+  qint32 xMin = -999999999, yMin = -999999999, zMin = -999999999;
+  qint32 xMax = 999999999, yMax = 999999999, zMax = 999999999;
+  QVector<qint32> data( count * 3 );
+  for ( int i = 0; i < count; ++i )
+  {
+    QByteArray bytes = f.read( stride );
+    // WHY??? X,Y,Z are int32 values stored as doubles
+    double *bytesD = (double*) bytes.constData();
+    data[i*3+0] = (bytesD[0]);
+    data[i*3+1] = (bytesD[1]);
+    data[i*3+2] = (bytesD[2]);
+
+    xMin = std::min( xMin, data[i*3+0]);
+    xMax = std::max( xMax, data[i*3+0]);
+    yMin = std::min( yMin, data[i*3+1]);
+    yMax = std::max( yMax, data[i*3+1]);
+    zMin = std::min( zMin, data[i*3+2]);
+    zMax = std::max( zMax, data[i*3+2]);
+  }
+  db = QgsPointCloudDataBounds(xMin, xMax, yMin, yMax, zMin, zMax);
+  return data;
+}
+
+/* *************************************************************************************** */
+
+QByteArray decompressZtdStream( const QByteArray &dataCompressed )
+{
+  // NOTE: this is very primitive implementation because we expect the uncompressed
+  // data will be always less than 10 MB
+
+  const int MAXSIZE=10000000;
+  QByteArray dataUncompressed;
+  dataUncompressed.resize( MAXSIZE );
+
+  ZSTD_DStream *strm = ZSTD_createDStream();
+  ZSTD_initDStream(strm);
+
+  ZSTD_inBuffer m_inBuf;
+  m_inBuf.src = reinterpret_cast<const void *>(dataCompressed.constData());
+  m_inBuf.size = dataCompressed.size();
+  m_inBuf.pos = 0;
+
+  ZSTD_outBuffer outBuf { reinterpret_cast<void *>(dataUncompressed.data()), MAXSIZE, 0 };
+  size_t ret = ZSTD_decompressStream(strm, &outBuf, &m_inBuf);
+  Q_ASSERT (!ZSTD_isError(ret));
+  Q_ASSERT( outBuf.pos );
+  Q_ASSERT( outBuf.pos < outBuf.size );
+
+  ZSTD_freeDStream(strm);
+  dataUncompressed.resize(outBuf.pos);
+  return dataUncompressed;
+}
+
+QVector<qint32> QgsPointCloudDecoder::decompressZStandard(const QString& filename, QgsPointCloudDataBounds &db)
+{
+  Q_ASSERT( QFile::exists( filename ) );
+
+  QFile f( filename );
+  bool r = f.open(QIODevice::ReadOnly);
+  Q_ASSERT(r);
+
+  QByteArray dataCompressed = f.readAll();
+  QByteArray dataUncompressed = decompressZtdStream( dataCompressed );
+
+  // from here it's the same as "binary"
+
+  // WHY??? per-record should be 18 based on schema, not 46
+  int stride = 46; //18;
+  int count = dataUncompressed.size() / stride;
+  qint32 xMin = -999999999, yMin = -999999999, zMin = -999999999;
+  qint32 xMax = 999999999, yMax = 999999999, zMax = 999999999;
+
+  QVector<qint32> data( count * 3 );
+  const char *ptr = dataUncompressed.constData();
+  for ( int i = 0; i < count; ++i )
+  {
+    // WHY??? X,Y,Z are int32 values stored as doubles
+    double *bytesD = (double*) (ptr+stride*i);
+    data[i*3+0] = (bytesD[0]);
+    data[i*3+1] = (bytesD[1]);
+    data[i*3+2] = (bytesD[2]);
+
+    xMin = std::min( xMin, data[i*3+0]);
+    xMax = std::max( xMax, data[i*3+0]);
+    yMin = std::min( yMin, data[i*3+1]);
+    yMax = std::max( yMax, data[i*3+1]);
+    zMin = std::min( zMin, data[i*3+2]);
+    zMax = std::max( zMax, data[i*3+2]);
+  }
+  db = QgsPointCloudDataBounds(xMin, xMax, yMin, yMax, zMin, zMax);
+  return data;
+}
+
+/* *************************************************************************************** */
+
+QVector<qint32> QgsPointCloudDecoder::decompressLaz(const QString& filename, QgsPointCloudDataBounds &db)
+{
+  std::ifstream file(filename.toLatin1().constData(), std::ios::binary);
+  Q_ASSERT (file.good());
+
+  auto start = common::tick();
+
+  laszip::io::reader::file f(file);
+
+  size_t count = f.get_header().point_count;
+  char buf[256]; // a buffer large enough to hold our point
+
+  qint32 xMin = -999999999, yMin = -999999999, zMin = -999999999;
+  qint32 xMax = 999999999, yMax = 999999999, zMax = 999999999;
+  QVector<qint32> data( count * 3 );
+
+  for(size_t i = 0 ; i < count ; i ++) {
+          f.readPoint(buf); // read the point out
+          laszip::formats::las::point10 p = laszip::formats::packers<laszip::formats::las::point10>::unpack(buf);
+
+          data[i*3+0] = p.x ;
+          data[i*3+1] = p.y ;
+          data[i*3+2] = p.z ;
+
+          xMin = std::min( xMin, data[i*3+0]);
+          xMax = std::max( xMax, data[i*3+0]);
+          yMin = std::min( yMin, data[i*3+1]);
+          yMax = std::max( yMax, data[i*3+1]);
+          zMin = std::min( zMin, data[i*3+2]);
+          zMax = std::max( zMax, data[i*3+2]);
+  }
+  db = QgsPointCloudDataBounds(xMin, xMax, yMin, yMax, zMin, zMax);
+  float t = common::since(start);
+  std::cout << "LAZ-PERF Read through the points in " << t << " seconds." << std::endl;
+}
diff --git a/src/core/pointcloud/qgspointclouddecoder.h b/src/core/pointcloud/qgspointclouddecoder.h
@@ -0,0 +1,41 @@
+/***************************************************************************
+                         qgspointclouddecoder.h
+                         --------------------
+    begin                : October 2020
+    copyright            : (C) 2020 by Peter Petrik
+    email                : zilolv 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.                                   *
+ *                                                                         *
+ ***************************************************************************/
+
+#ifndef QGSPOINTCLOUDDECODER_H
+#define QGSPOINTCLOUDDECODER_H
+
+
+#include "qgis_core.h"
+#include "qgis_sip.h"
+
+#define SIP_NO_FILE
+
+#include <QVector>
+#include <QString>
+
+class QgsPointCloudDataBounds;
+class QgsVector3D;
+
+namespace QgsPointCloudDecoder
+{
+  QVector<qint32> decompressBinary(const QString& filename, QgsPointCloudDataBounds &db);
+  QVector<qint32> decompressZStandard(const QString& filename, QgsPointCloudDataBounds &db);
+  QVector<qint32> decompressLaz(const QString& filename, QgsPointCloudDataBounds &db);
+};
+
+
+#endif // QGSPOINTCLOUDDECODER_H