Skip to content

Commit

Permalink
Make coordinates in point cloud chunks independent from scene origin
Browse files Browse the repository at this point in the history 
Until now, coordinates of chunks of point clouds were using "world" coordinates,
i.e. they were relative to the 3D scene origin point. As a part of preparation
for large scene support (QEP #301), we need coordinates inside chunks to be
independent from the scene origin point. For each point cloud chunk, we pick
origin point as a corner of the chunk's 3D box, which should guarantee that
coordinates within chunk (floats) will have reasonable numerical precision.

I have also used the opportunity to simplify the code that handles point cloud
chunk coordinates: rather than having flipped coordinates (x,-z,y) inside chunks,
we are using map coordinates without flipping (x,y,z), only relative to the origin
point of the chunk. The flipping of axes is handled by QTransform attached
to each chunk. (Hopefully when all chunked entity implementations get updated,
we will be able to remove all the axis flipping everywhere by just removing
the rotation in associated QTransforms.)

There's also a fix for a crash, when point cloud uses coloring based on classification
and triangulation is enabled (the pointSize array was not getting populated for parent
nodes).
  • Loading branch information
@wonder-sk
wonder-sk committed Oct 2, 2024
1 parent 7d20df7 commit bc40b15
Show file tree
Hide file tree
Showing 2 changed files with 125 additions and 60 deletions.
168 changes: 115 additions & 53 deletions src/3d/symbols/qgspointcloud3dsymbol_p.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -54,6 +54,28 @@ typedef Qt3DCore::QGeometry Qt3DQGeometry;

#include <delaunator.hpp>

// pick a point that we'll use as origin for coordinates for this node's points
static QgsVector3D originFromNodeBounds( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context, const QgsPointCloudBlock *block )
{
const QgsVector3D blockScale = block->scale();
const QgsVector3D blockOffset = block->offset();

QgsPointCloudDataBounds bounds = pc->nodeBounds( n );
double nodeOriginX = bounds.xMin() * blockScale.x() + blockOffset.x();
double nodeOriginY = bounds.yMin() * blockScale.y() + blockOffset.y();
double nodeOriginZ = ( bounds.zMin() * blockScale.z() + blockOffset.z() ) * context.zValueScale() + context.zValueFixedOffset();
try
{
context.coordinateTransform().transformInPlace( nodeOriginX, nodeOriginY, nodeOriginZ );
}
catch ( QgsCsException & )
{
QgsDebugError( QStringLiteral( "Error transforming node origin point" ) );
}
return QgsVector3D( nodeOriginX, nodeOriginY, nodeOriginZ );
}


QgsPointCloud3DGeometry::QgsPointCloud3DGeometry( Qt3DCore::QNode *parent, const QgsPointCloud3DSymbolHandler::PointData &data, unsigned int byteStride )
: Qt3DQGeometry( parent )
, mPositionAttribute( new Qt3DQAttribute( this ) )
Expand Down Expand Up @@ -241,6 +263,7 @@ void QgsClassificationPointCloud3DGeometry::makeVertexBuffer( const QgsPointClou
float *rawVertexArray = reinterpret_cast<float *>( vertexBufferData.data() );
int idx = 0;
Q_ASSERT( data.positions.size() == data.parameter.size() );
Q_ASSERT( data.positions.size() == data.pointSizes.size() );
for ( int i = 0; i < data.positions.size(); ++i )
{
rawVertexArray[idx++] = data.positions.at( i ).x();
Expand Down Expand Up @@ -281,8 +304,12 @@ void QgsPointCloud3DSymbolHandler::makeEntity( Qt3DCore::QEntity *parent, const
}
gr->setGeometry( geom );

// Transform
// Transform: chunks are using coordinates relative to chunk origin, with X,Y,Z axes being the same
// as map coordinates, so we need to rotate and translate entities to get them into world coordinates
Qt3DCore::QTransform *tr = new Qt3DCore::QTransform;
QVector3D nodeTranslation = ( out.positionsOrigin - context.origin() ).toVector3D();
tr->setRotation( QQuaternion::fromAxisAndAngle( QVector3D( 1, 0, 0 ), -90 ) ); // flip map (x,y,z) to world (x,z,-y)
tr->setTranslation( QVector3D( nodeTranslation.x(), nodeTranslation.z(), -nodeTranslation.y() ) );

// Material
Qt3DRender::QMaterial *mat = new Qt3DRender::QMaterial;
Expand Down Expand Up @@ -337,63 +364,64 @@ std::vector<double> QgsPointCloud3DSymbolHandler::getVertices( QgsPointCloudInde

bool hasColorData = !outNormal.colors.empty();
bool hasParameterData = !outNormal.parameter.empty();
bool hasPointSizeData = !outNormal.pointSizes.empty();

// first, get the points of the concerned node
std::vector<double> vertices( outNormal.positions.size() * 2 );
size_t idx = 0;
for ( int i = 0; i < outNormal.positions.size(); ++i )
{
vertices[idx++] = outNormal.positions.at( i ).x();
vertices[idx++] = outNormal.positions.at( i ).z();
vertices[idx++] = -outNormal.positions.at( i ).y(); // flipping y to have correctly oriented triangles from delaunator
}

// next, we also need all points of all parents nodes to make the triangulation (also external points)
IndexedPointCloudNode parentNode = n.parentNode();

int properPointsCount = outNormal.positions.count();
while ( parentNode.d() >= 0 )
{
processNode( pc, parentNode, context );
parentNode = parentNode.parentNode();
}

PointData filteredExtraPointData;

double span = pc->span();
//factor to take account of the density of the point to calculate extension of the bounding box
// with a usual value span = 128, bounding box is extended by 12.5 % on each side.
double extraBoxFactor = 16 / span;
double extraX = extraBoxFactor * bbox.xExtent();
double extraZ = extraBoxFactor * bbox.zExtent();
double extraY = extraBoxFactor * bbox.yExtent();

// We keep all points in vertical direction to avoid odd triangulation if points are isolated on top
const QgsAABB extendedBBox( bbox.xMin - extraX, -std::numeric_limits<float>::max(), bbox.zMin - extraZ, bbox.xMax + extraX, std::numeric_limits<float>::max(), bbox.zMax + extraZ );
const QgsAABB extendedBBox( bbox.xMin - extraX, bbox.yMin - extraY, -std::numeric_limits<float>::max(), bbox.xMax + extraX, bbox.yMax + extraY, std::numeric_limits<float>::max() );

for ( int i = properPointsCount; i < outNormal.positions.count(); ++i )
PointData filteredExtraPointData;
while ( parentNode.d() >= 0 )
{
const QVector3D pos = outNormal.positions.at( i );
if ( extendedBBox.intersects( pos.x(), pos.y(), pos.z() ) )
PointData outputParent;
processNode( pc, parentNode, context, &outputParent );

// the "main" chunk and each parent chunks have their origins
QVector3D originDifference = ( outputParent.positionsOrigin - outNormal.positionsOrigin ).toVector3D();

for ( int i = 0; i < outputParent.positions.count(); ++i )
{
filteredExtraPointData.positions.append( pos );
vertices.push_back( pos.x() );
vertices.push_back( pos.z() );

if ( hasColorData )
filteredExtraPointData.colors.append( outNormal.colors.at( i ) );
if ( hasParameterData )
filteredExtraPointData.parameter.append( outNormal.parameter.at( i ) );
const QVector3D pos = outputParent.positions.at( i ) + originDifference;
if ( extendedBBox.intersects( pos.x(), pos.y(), pos.z() ) )
{
filteredExtraPointData.positions.append( pos );
vertices.push_back( pos.x() );
vertices.push_back( -pos.y() ); // flipping y to have correctly oriented triangles from delaunator

if ( hasColorData )
filteredExtraPointData.colors.append( outputParent.colors.at( i ) );
if ( hasParameterData )
filteredExtraPointData.parameter.append( outputParent.parameter.at( i ) );
if ( hasPointSizeData )
filteredExtraPointData.pointSizes.append( outputParent.pointSizes.at( i ) );
}
}
}

outNormal.positions.resize( properPointsCount );
if ( hasColorData )
outNormal.colors.resize( properPointsCount );
if ( hasParameterData )
outNormal.parameter.resize( properPointsCount );
parentNode = parentNode.parentNode();
}

outNormal.positions.append( filteredExtraPointData.positions );
outNormal.colors.append( filteredExtraPointData.colors );
outNormal.parameter.append( filteredExtraPointData.parameter );
outNormal.pointSizes.append( filteredExtraPointData.pointSizes );

return vertices;
}
Expand Down Expand Up @@ -456,13 +484,13 @@ void QgsPointCloud3DSymbolHandler::filterTriangles( const std::vector<size_t> &t
const QVector3D pos2 = outNormal.positions.at( triangleIndexes.at( i + ( j + 1 ) % 3 ) );

if ( verticalFilter )
verticalSkip |= std::fabs( pos.y() - pos2.y() ) > verticalThreshold;
verticalSkip |= std::fabs( pos.z() - pos2.z() ) > verticalThreshold;

if ( horizontalFilter && ! verticalSkip )
{
// filter only in the horizontal plan, it is a 2.5D triangulation.
horizontalSkip |= sqrt( std::pow( pos.x() - pos2.x(), 2 ) +
std::pow( pos.z() - pos2.z(), 2 ) ) > horizontalThreshold;
std::pow( pos.y() - pos2.y(), 2 ) ) > horizontalThreshold;
}

if ( horizontalSkip || verticalSkip )
Expand Down Expand Up @@ -496,11 +524,25 @@ void QgsPointCloud3DSymbolHandler::triangulate( QgsPointCloudIndex *pc, const In
if ( outNormal.positions.isEmpty() )
return;

// The bbox we get is in world coordinates, we need to transform it to map coordinates
// (flip axes and add scene origin vector), but relative to the chunk's origin (subtract it)
// because that's the coordinate system used within the chunk
QgsBox3D boxRelativeToChunkOrigin(
bbox.xMin + context.origin().x() - outNormal.positionsOrigin.x(),
-bbox.zMax + context.origin().y() - outNormal.positionsOrigin.y(),
bbox.yMin + context.origin().z() - outNormal.positionsOrigin.z(),
bbox.xMax + context.origin().x() - outNormal.positionsOrigin.x(),
-bbox.zMin + context.origin().y() - outNormal.positionsOrigin.y(),
bbox.yMax + context.origin().z() - outNormal.positionsOrigin.z() );
QgsAABB aabbRelativeToChunkOrigin(
boxRelativeToChunkOrigin.xMinimum(), boxRelativeToChunkOrigin.yMinimum(), boxRelativeToChunkOrigin.zMinimum(),
boxRelativeToChunkOrigin.xMaximum(), boxRelativeToChunkOrigin.yMaximum(), boxRelativeToChunkOrigin.zMaximum() );

// Triangulation happens here
std::unique_ptr<delaunator::Delaunator> triangulation;
try
{
std::vector<double> vertices = getVertices( pc, n, context, bbox );
std::vector<double> vertices = getVertices( pc, n, context, aabbRelativeToChunkOrigin );
triangulation.reset( new delaunator::Delaunator( vertices ) );
}
catch ( std::exception &e )
Expand All @@ -515,7 +557,7 @@ void QgsPointCloud3DSymbolHandler::triangulate( QgsPointCloudIndex *pc, const In
const std::vector<size_t> &triangleIndexes = triangulation->triangles;

calculateNormals( triangleIndexes );
filterTriangles( triangleIndexes, context, bbox );
filterTriangles( triangleIndexes, context, aabbRelativeToChunkOrigin );
}

std::unique_ptr<QgsPointCloudBlock> QgsPointCloud3DSymbolHandler::pointCloudBlock( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloudRequest &request, const QgsPointCloud3DRenderContext &context )
Expand Down Expand Up @@ -561,7 +603,7 @@ bool QgsSingleColorPointCloud3DSymbolHandler::prepare( const QgsPointCloud3DRend
return true;
}

void QgsSingleColorPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context )
void QgsSingleColorPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context, PointData *output )
{
QgsPointCloudAttributeCollection attributes;
attributes.push_back( QgsPointCloudAttribute( QStringLiteral( "X" ), QgsPointCloudAttribute::Int32 ) );
Expand All @@ -585,6 +627,11 @@ void QgsSingleColorPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *p
const QgsCoordinateTransform coordinateTransform = context.coordinateTransform();
bool alreadyPrintedDebug = false;

if ( !output )
output = &outNormal;

output->positionsOrigin = originFromNodeBounds( pc, n, context, block.get() );

for ( int i = 0; i < count; ++i )
{
if ( context.isCanceled() )
Expand All @@ -609,8 +656,8 @@ void QgsSingleColorPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *p
alreadyPrintedDebug = true;
}
}
const QgsVector3D point = context.mapToWorldCoordinates( QgsVector3D( x, y, z ) );
outNormal.positions.push_back( QVector3D( static_cast<float>( point.x() ), static_cast<float>( point.y() ), static_cast<float>( point.z() ) ) );
const QgsVector3D point = QgsVector3D( x, y, z ) - output->positionsOrigin;
output->positions.push_back( point.toVector3D() );
}
}

Expand All @@ -636,7 +683,7 @@ bool QgsColorRampPointCloud3DSymbolHandler::prepare( const QgsPointCloud3DRender
return true;
}

void QgsColorRampPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context )
void QgsColorRampPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context, PointData *output )
{
QgsPointCloudAttributeCollection attributes;
const int xOffset = 0;
Expand Down Expand Up @@ -705,6 +752,11 @@ void QgsColorRampPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc,
const QgsVector3D blockScale = block->scale();
const QgsVector3D blockOffset = block->offset();

if ( !output )
output = &outNormal;

output->positionsOrigin = originFromNodeBounds( pc, n, context, block.get() );

for ( int i = 0; i < count; ++i )
{
if ( context.isCanceled() )
Expand All @@ -729,20 +781,20 @@ void QgsColorRampPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc,
alreadyPrintedDebug = true;
}
}
const QgsVector3D point = context.mapToWorldCoordinates( QgsVector3D( x, y, z ) );
outNormal.positions.push_back( QVector3D( static_cast<float>( point.x() ), static_cast<float>( point.y() ), static_cast<float>( point.z() ) ) );
const QgsVector3D point = QgsVector3D( x, y, z ) - output->positionsOrigin;
output->positions.push_back( point.toVector3D() );

if ( attrIsX )
outNormal.parameter.push_back( x );
output->parameter.push_back( x );
else if ( attrIsY )
outNormal.parameter.push_back( y );
output->parameter.push_back( y );
else if ( attrIsZ )
outNormal.parameter.push_back( z );
output->parameter.push_back( z );
else
{
float iParam = 0.0f;
context.getAttribute( ptr, i * recordSize + attributeOffset, attributeType, iParam );
outNormal.parameter.push_back( iParam );
output->parameter.push_back( iParam );
}
}
}
Expand All @@ -769,7 +821,7 @@ bool QgsRGBPointCloud3DSymbolHandler::prepare( const QgsPointCloud3DRenderContex
return true;
}

void QgsRGBPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context )
void QgsRGBPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context, PointData *output )
{
QgsPointCloudAttributeCollection attributes;
attributes.push_back( QgsPointCloudAttribute( QStringLiteral( "X" ), QgsPointCloudAttribute::Int32 ) );
Expand Down Expand Up @@ -822,6 +874,11 @@ void QgsRGBPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const
const bool useBlueContrastEnhancement = blueContrastEnhancement && blueContrastEnhancement->contrastEnhancementAlgorithm() != QgsContrastEnhancement::NoEnhancement;
const bool useGreenContrastEnhancement = greenContrastEnhancement && greenContrastEnhancement->contrastEnhancementAlgorithm() != QgsContrastEnhancement::NoEnhancement;

if ( !output )
output = &outNormal;

output->positionsOrigin = originFromNodeBounds( pc, n, context, block.get() );

int ir = 0;
int ig = 0;
int ib = 0;
Expand All @@ -848,7 +905,7 @@ void QgsRGBPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const
alreadyPrintedDebug = true;
}
}
const QgsVector3D point = context.mapToWorldCoordinates( QgsVector3D( x, y, z ) );
const QgsVector3D point = QgsVector3D( x, y, z ) - output->positionsOrigin;

QVector3D color( 0.0f, 0.0f, 0.0f );

Expand Down Expand Up @@ -882,8 +939,8 @@ void QgsRGBPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const
color.setY( ig / 255.0f );
color.setZ( ib / 255.0f );

outNormal.positions.push_back( QVector3D( static_cast<float>( point.x() ), static_cast<float>( point.y() ), static_cast<float>( point.z() ) ) );
outNormal.colors.push_back( color );
output->positions.push_back( point.toVector3D() );
output->colors.push_back( color );
}
}

Expand All @@ -909,7 +966,7 @@ bool QgsClassificationPointCloud3DSymbolHandler::prepare( const QgsPointCloud3DR
return true;
}

void QgsClassificationPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context )
void QgsClassificationPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex *pc, const IndexedPointCloudNode &n, const QgsPointCloud3DRenderContext &context, PointData *output )
{
QgsPointCloudAttributeCollection attributes;
const int xOffset = 0;
Expand Down Expand Up @@ -987,6 +1044,11 @@ void QgsClassificationPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex
context.symbol() ? context.symbol()->pointSize() : 1.0 );
}

if ( !output )
output = &outNormal;

output->positionsOrigin = originFromNodeBounds( pc, n, context, block.get() );

const QSet<int> filteredOutValues = context.getFilteredOutValues();
for ( int i = 0; i < count; ++i )
{
Expand All @@ -1012,7 +1074,7 @@ void QgsClassificationPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex
alreadyPrintedDebug = true;
}
}
const QgsVector3D point = context.mapToWorldCoordinates( QgsVector3D( x, y, z ) );
const QgsVector3D point = QgsVector3D( x, y, z ) - output->positionsOrigin;
float iParam = 0.0f;
if ( attrIsX )
iParam = x;
Expand All @@ -1026,12 +1088,12 @@ void QgsClassificationPointCloud3DSymbolHandler::processNode( QgsPointCloudIndex
if ( filteredOutValues.contains( ( int ) iParam ) ||
! categoriesValues.contains( ( int ) iParam ) )
continue;
outNormal.positions.push_back( QVector3D( static_cast<float>( point.x() ), static_cast<float>( point.y() ), static_cast<float>( point.z() ) ) );
output->positions.push_back( point.toVector3D() );

// find iParam actual index in the categories list
float iParam2 = categoriesValues.indexOf( ( int )iParam ) + 1;
outNormal.parameter.push_back( iParam2 );
outNormal.pointSizes.push_back( categoriesPointSizes.value( ( int ) iParam ) );
output->parameter.push_back( iParam2 );
output->pointSizes.push_back( categoriesPointSizes.value( ( int ) iParam ) );
}
}

Expand Down
Loading

0 comments on commit bc40b15

Please sign in to comment.