Skip to content

Commit

Permalink
read CSR and bounding box
Browse files Browse the repository at this point in the history
  • Loading branch information
@PeterPetrik @nyalldawson
PeterPetrik authored and nyalldawson committed Oct 26, 2020
1 parent 8a42c57 commit 1d4b8e7
Show file tree
Hide file tree
Showing 7 changed files with 119 additions and 138 deletions.
98 changes: 35 additions & 63 deletions src/core/pointcloud/qgspointclouddecoder.cpp
View file
Expand Up @@ -26,37 +26,27 @@
#include "laz-perf/io.hpp"
#include "laz-perf/common/common.hpp"

QVector<qint32> QgsPointCloudDecoder::decompressBinary(const QString& filename, QgsPointCloudDataBounds &db)
QVector<qint32> QgsPointCloudDecoder::decompressBinary( const QString &filename )
{
Q_ASSERT( QFile::exists( filename ) );

QFile f( filename );
bool r = f.open(QIODevice::ReadOnly);
Q_ASSERT(r);
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]);
double *bytesD = ( double * ) bytes.constData();
data[i * 3 + 0] = ( bytesD[0] );
data[i * 3 + 1] = ( bytesD[1] );
data[i * 3 + 2] = ( bytesD[2] );
}
db = QgsPointCloudDataBounds(xMin, xMax, yMin, yMax, zMin, zMax);
return data;
}

Expand All @@ -67,36 +57,36 @@ 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;
const int MAXSIZE = 10000000;
QByteArray dataUncompressed;
dataUncompressed.resize( MAXSIZE );

ZSTD_DStream *strm = ZSTD_createDStream();
ZSTD_initDStream(strm);
ZSTD_initDStream( strm );

ZSTD_inBuffer m_inBuf;
m_inBuf.src = reinterpret_cast<const void *>(dataCompressed.constData());
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));
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);
ZSTD_freeDStream( strm );
dataUncompressed.resize( outBuf.pos );
return dataUncompressed;
}

QVector<qint32> QgsPointCloudDecoder::decompressZStandard(const QString& filename, QgsPointCloudDataBounds &db)
QVector<qint32> QgsPointCloudDecoder::decompressZStandard( const QString &filename )
{
Q_ASSERT( QFile::exists( filename ) );

QFile f( filename );
bool r = f.open(QIODevice::ReadOnly);
Q_ASSERT(r);
bool r = f.open( QIODevice::ReadOnly );
Q_ASSERT( r );

QByteArray dataCompressed = f.readAll();
QByteArray dataUncompressed = decompressZtdStream( dataCompressed );
Expand All @@ -106,64 +96,46 @@ QVector<qint32> QgsPointCloudDecoder::decompressZStandard(const QString& filenam
// 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]);
double *bytesD = ( double * )( ptr + stride * i );
data[i * 3 + 0] = ( bytesD[0] );
data[i * 3 + 1] = ( bytesD[1] );
data[i * 3 + 2] = ( bytesD[2] );
}
db = QgsPointCloudDataBounds(xMin, xMax, yMin, yMax, zMin, zMax);
return data;
}

/* *************************************************************************************** */

QVector<qint32> QgsPointCloudDecoder::decompressLaz(const QString& filename, QgsPointCloudDataBounds &db)
QVector<qint32> QgsPointCloudDecoder::decompressLaz( const QString &filename )
{
std::ifstream file(filename.toLatin1().constData(), std::ios::binary);
Q_ASSERT (file.good());
std::ifstream file( filename.toLatin1().constData(), std::ios::binary );
Q_ASSERT( file.good() );

auto start = common::tick();

laszip::io::reader::file f(file);
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 ;
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 );

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]);
data[i * 3 + 0] = p.x ;
data[i * 3 + 1] = p.y ;
data[i * 3 + 2] = p.z ;
}
db = QgsPointCloudDataBounds(xMin, xMax, yMin, yMax, zMin, zMax);
float t = common::since(start);
float t = common::since( start );
std::cout << "LAZ-PERF Read through the points in " << t << " seconds." << std::endl;
}
9 changes: 3 additions & 6 deletions src/core/pointcloud/qgspointclouddecoder.h
View file
Expand Up @@ -27,14 +27,11 @@
#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);
QVector<qint32> decompressBinary( const QString &filename );
QVector<qint32> decompressZStandard( const QString &filename );
QVector<qint32> decompressLaz( const QString &filename );
};


Expand Down
37 changes: 31 additions & 6 deletions src/core/pointcloud/qgspointcloudindex.cpp
View file
Expand Up @@ -116,6 +116,7 @@ bool QgsPointCloudIndex::load( const QString &fileName )
return false;

const QDir directory = QFileInfo( fileName ).absoluteDir();
mDirectory = directory.absolutePath();

QByteArray dataJson = f.readAll();
QJsonParseError err;
Expand All @@ -133,10 +134,19 @@ bool QgsPointCloudIndex::load( const QString &fileName )

mSpan = doc["span"].toInt();

// WKT
QJsonObject srs = doc["srs"].toObject();
mWkt = srs["wkt"].toString();

// rectangular
QJsonArray bounds = doc["bounds"].toArray();
if ( bounds.size() != 6 )
return false;

QJsonArray bounds_conforming = doc["boundsConforming"].toArray();
if ( bounds.size() != 6 )
return false;

QJsonArray schemaArray = doc["schema"].toArray();

for ( QJsonValue schemaItem : schemaArray )
Expand Down Expand Up @@ -201,7 +211,7 @@ bool QgsPointCloudIndex::load( const QString &fileName )

// load hierarchy

QFile fH( directory.filePath( "/ept-hierarchy/0-0-0-0.json" ) );
QFile fH( QStringLiteral( "%1/ept-hierarchy/0-0-0-0.json" ).arg( mDirectory ) );
if ( !fH.open( QIODevice::ReadOnly ) )
return false;

Expand Down Expand Up @@ -242,28 +252,28 @@ QList<IndexedPointCloudNode> QgsPointCloudIndex::children( const IndexedPointClo
return lst;
}

QVector<qint32> QgsPointCloudIndex::nodePositionDataAsInt32( const IndexedPointCloudNode &n, QgsVector3D, QgsVector3D, QgsPointCloudDataBounds &db )
QVector<qint32> QgsPointCloudIndex::nodePositionDataAsInt32( const IndexedPointCloudNode &n )
{
Q_ASSERT( mHierarchy.contains( n ) );
int count = mHierarchy[n];
// int count = mHierarchy[n];

if ( mDataType == "binary" )
{
QString filename = QString( "%1/ept-data/%2.bin" ).arg( mDirectory ).arg( n.toString() );
Q_ASSERT( QFile::exists( filename ) );
return QgsPointCloudDecoder::decompressBinary( filename, db );
return QgsPointCloudDecoder::decompressBinary( filename );
}
else if ( mDataType == "zstandard" )
{
QString filename = QString( "%1/ept-data/%2.zst" ).arg( mDirectory ).arg( n.toString() );
Q_ASSERT( QFile::exists( filename ) );
return QgsPointCloudDecoder::decompressBinary( filename, db );
return QgsPointCloudDecoder::decompressBinary( filename );
}
else // if ( mDataType == "laz" )
{
QString filename = QString( "%1/ept-data/%2.laz" ).arg( mDirectory ).arg( n.toString() );
Q_ASSERT( QFile::exists( filename ) );
return QgsPointCloudDecoder::decompressBinary( filename, db );
return QgsPointCloudDecoder::decompressBinary( filename );
}
}

Expand All @@ -285,3 +295,18 @@ QgsPointCloudDataBounds QgsPointCloudIndex::nodeBounds( const IndexedPointCloudN
QgsPointCloudDataBounds db( xMin, xMax, yMin, yMax, zMin, zMax );
return db;
}

QString QgsPointCloudIndex::wkt() const
{
return mWkt;
}

QgsVector3D QgsPointCloudIndex::scale() const
{
return mScale;
}

QgsVector3D QgsPointCloudIndex::offset() const
{
return mOffset;
}
9 changes: 8 additions & 1 deletion src/core/pointcloud/qgspointcloudindex.h
View file
Expand Up @@ -107,9 +107,14 @@ class CORE_EXPORT QgsPointCloudIndex: public QObject

QList<IndexedPointCloudNode> children( const IndexedPointCloudNode &n );

QVector<qint32> nodePositionDataAsInt32( const IndexedPointCloudNode &n, QgsVector3D scale, QgsVector3D offset, QgsPointCloudDataBounds &db );
QVector<qint32> nodePositionDataAsInt32( const IndexedPointCloudNode &n );

QgsPointCloudDataBounds nodeBounds( const IndexedPointCloudNode &n );
QString wkt() const;

QgsVector3D scale() const;

QgsVector3D offset() const;

private:
QString mDirectory;
Expand All @@ -120,6 +125,8 @@ class CORE_EXPORT QgsPointCloudIndex: public QObject
QgsVector3D mOffset; //!< Offset of our int32 coordinates compared to CRS coords
QgsPointCloudDataBounds mRootBounds; //!< Bounds of the root node's cube (in int32 coordinates)
int mSpan; //!< Number of points in one direction in a single node

QString mWkt;
};


Expand Down
23 changes: 20 additions & 3 deletions src/core/pointcloud/qgspointcloudlayer.cpp
View file
Expand Up @@ -21,7 +21,7 @@
#include "qgsrectangle.h"

QgsPointCloudLayer::QgsPointCloudLayer( const QString &path, const QString &baseName )
: QgsMapLayer( QgsMapLayerType::PointCloudLayer, path, baseName )
: QgsMapLayer( QgsMapLayerType::PointCloudLayer, baseName, path )
{
setValid( loadDataSource() );
}
Expand All @@ -37,7 +37,19 @@ QgsPointCloudLayer *QgsPointCloudLayer::clone() const

QgsRectangle QgsPointCloudLayer::extent() const
{
return QgsRectangle();
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 );
}

QgsMapLayerRenderer *QgsPointCloudLayer::createMapRenderer( QgsRenderContext &rendererContext )
Expand Down Expand Up @@ -101,5 +113,10 @@ QgsPointCloudIndex *QgsPointCloudLayer::pointCloudIndex() const
bool QgsPointCloudLayer::loadDataSource()
{
mPointCloudIndex.reset( new QgsPointCloudIndex() );
return mPointCloudIndex->load( source() );
bool success = mPointCloudIndex->load( source() );

if ( success )
setCrs( QgsCoordinateReferenceSystem::fromWkt( mPointCloudIndex->wkt() ) );

return success;
}

0 comments on commit 1d4b8e7

Please sign in to comment.