Skip to content

Commit

Permalink
Fix elevation profile on virtual point cloud (qgis#59073)
Browse files Browse the repository at this point in the history 
* Fix elevation profile on virtual point cloud

(fix qgis#54728)

* Update src/core/pointcloud/qgspointcloudlayerprofilegenerator.cpp

---------

Co-authored-by: Nyall Dawson <nyall.dawson@gmail.com>
  • Loading branch information
@dvdkon @nyalldawson
dvdkon and nyalldawson authored Oct 24, 2024
1 parent 400cf13 commit de9377a
Showing 1 changed file with 81 additions and 49 deletions.
130 changes: 81 additions & 49 deletions src/core/pointcloud/qgspointcloudlayerprofilegenerator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -389,12 +389,23 @@ bool QgsPointCloudLayerProfileGenerator::generateProfile( const QgsProfileGenera
// this is not AT ALL thread safe, but it's what QgsPointCloudLayerRenderer does !
// TODO: fix when QgsPointCloudLayerRenderer is made thread safe to use same approach

QgsPointCloudIndex *pc = mLayer->dataProvider()->index();
if ( !pc || !pc->isValid() )
QVector<QgsPointCloudIndex *> indexes;
QgsPointCloudIndex *mainIndex = mLayer->dataProvider()->index();
if ( mainIndex && mainIndex->isValid() )
indexes.append( mainIndex );

// Gather all relevant sub-indexes
const QgsRectangle profileCurveBbox = mProfileCurve->boundingBox();
for ( const QgsPointCloudSubIndex &subidx : mLayer->dataProvider()->subIndexes() )
{
return false;
QgsPointCloudIndex *index = subidx.index();
if ( index && index->isValid() && subidx.polygonBounds().intersects( profileCurveBbox ) )
indexes.append( subidx.index() );
}

if ( indexes.empty() )
return false;

const double startDistanceOffset = std::max( !context.distanceRange().isInfinite() ? context.distanceRange().lower() : 0, 0.0 );
const double endDistance = context.distanceRange().upper();

Expand Down Expand Up @@ -433,9 +444,13 @@ bool QgsPointCloudLayerProfileGenerator::generateProfile( const QgsProfileGenera
mSearchGeometryInLayerCrsGeometryEngine->prepareGeometry();
mMaxSearchExtentInLayerCrs = mSearchGeometryInLayerCrs->boundingBox();

const IndexedPointCloudNode root = pc->root();

double maximumErrorPixels = context.convertDistanceToPixels( mMaximumScreenError, mMaximumScreenErrorUnit );
if ( maximumErrorPixels < 0.0 )
{
QgsDebugError( QStringLiteral( "Invalid maximum error in pixels" ) );
return false;
}

const double toleranceInPixels = context.convertDistanceToPixels( mTolerance, Qgis::RenderUnit::MapUnits );
// ensure that the maximum error is compatible with the tolerance size -- otherwise if the tolerance size
// is much smaller than the maximum error, we don't dig deep enough into the point cloud nodes to find
Expand All @@ -444,42 +459,6 @@ bool QgsPointCloudLayerProfileGenerator::generateProfile( const QgsProfileGenera
if ( toleranceInPixels / 4 < maximumErrorPixels )
maximumErrorPixels = toleranceInPixels / 4;

const QgsRectangle rootNodeExtentLayerCoords = pc->nodeMapExtent( root );
QgsRectangle rootNodeExtentInCurveCrs;
try
{
QgsCoordinateTransform extentTransform = mLayerToTargetTransform;
extentTransform.setBallparkTransformsAreAppropriate( true );
rootNodeExtentInCurveCrs = extentTransform.transformBoundingBox( rootNodeExtentLayerCoords );
}
catch ( QgsCsException & )
{
QgsDebugError( QStringLiteral( "Could not transform node extent to curve CRS" ) );
rootNodeExtentInCurveCrs = rootNodeExtentLayerCoords;
}

const double rootErrorInMapCoordinates = rootNodeExtentInCurveCrs.width() / pc->span(); // in curve coords

const double mapUnitsPerPixel = context.mapUnitsPerDistancePixel();
if ( ( rootErrorInMapCoordinates < 0.0 ) || ( mapUnitsPerPixel < 0.0 ) || ( maximumErrorPixels < 0.0 ) )
{
QgsDebugError( QStringLiteral( "invalid screen error" ) );
return false;
}
double rootErrorPixels = rootErrorInMapCoordinates / mapUnitsPerPixel; // in pixels
const QVector<IndexedPointCloudNode> nodes = traverseTree( pc, pc->root(), maximumErrorPixels, rootErrorPixels, context.elevationRange() );
if ( nodes.empty() )
{
return false;
}

const double rootErrorInLayerCoordinates = rootNodeExtentLayerCoords.width() / pc->span();
const double maxErrorInMapCoordinates = maximumErrorPixels * mapUnitsPerPixel;

mResults = std::make_unique< QgsPointCloudLayerProfileResults >();
mResults->copyPropertiesFromGenerator( this );
mResults->mMaxErrorInLayerCoordinates = maxErrorInMapCoordinates * rootErrorInLayerCoordinates / rootErrorInMapCoordinates;

QgsPointCloudRequest request;
QgsPointCloudAttributeCollection attributes;
attributes.push_back( QgsPointCloudAttribute( QStringLiteral( "X" ), QgsPointCloudAttribute::Int32 ) );
Expand Down Expand Up @@ -515,22 +494,75 @@ bool QgsPointCloudLayerProfileGenerator::generateProfile( const QgsProfileGenera

request.setAttributes( attributes );

switch ( pc->accessType() )
mResults = std::make_unique< QgsPointCloudLayerProfileResults >();
mResults->copyPropertiesFromGenerator( this );
mResults->mMaxErrorInLayerCoordinates = 0;

QgsCoordinateTransform extentTransform = mLayerToTargetTransform;
extentTransform.setBallparkTransformsAreAppropriate( true );

const double mapUnitsPerPixel = context.mapUnitsPerDistancePixel();
if ( mapUnitsPerPixel < 0.0 )
{
QgsDebugError( QStringLiteral( "Invalid map units per pixel ratio" ) );
return false;
}

for ( QgsPointCloudIndex *pc : std::as_const( indexes ) )
{
case QgsPointCloudIndex::AccessType::Local:
const IndexedPointCloudNode root = pc->root();
const QgsRectangle rootNodeExtentLayerCoords = pc->nodeMapExtent( root );
QgsRectangle rootNodeExtentInCurveCrs;
try
{
visitNodesSync( nodes, pc, request, context.elevationRange() );
break;
rootNodeExtentInCurveCrs = extentTransform.transformBoundingBox( rootNodeExtentLayerCoords );
}
case QgsPointCloudIndex::AccessType::Remote:
catch ( QgsCsException & )
{
visitNodesAsync( nodes, pc, request, context.elevationRange() );
break;
QgsDebugError( QStringLiteral( "Could not transform node extent to curve CRS" ) );
rootNodeExtentInCurveCrs = rootNodeExtentLayerCoords;
}

const double rootErrorInMapCoordinates = rootNodeExtentInCurveCrs.width() / pc->span(); // in curve coords
if ( rootErrorInMapCoordinates < 0.0 )
{
QgsDebugError( QStringLiteral( "Invalid root node error" ) );
return false;
}

double rootErrorPixels = rootErrorInMapCoordinates / mapUnitsPerPixel; // in pixels
const QVector<IndexedPointCloudNode> nodes = traverseTree( pc, pc->root(), maximumErrorPixels, rootErrorPixels, context.elevationRange() );

const double rootErrorInLayerCoordinates = rootNodeExtentLayerCoords.width() / pc->span();
const double maxErrorInMapCoordinates = maximumErrorPixels * mapUnitsPerPixel;

mResults->mMaxErrorInLayerCoordinates = std::max(
mResults->mMaxErrorInLayerCoordinates,
maxErrorInMapCoordinates * rootErrorInLayerCoordinates / rootErrorInMapCoordinates );

switch ( pc->accessType() )
{
case QgsPointCloudIndex::AccessType::Local:
{
visitNodesSync( nodes, pc, request, context.elevationRange() );
break;
}
case QgsPointCloudIndex::AccessType::Remote:
{
visitNodesAsync( nodes, pc, request, context.elevationRange() );
break;
}
}

if ( mFeedback->isCanceled() )
return false;
}

if ( mFeedback->isCanceled() )
if ( mGatheredPoints.empty() )
{
mResults = nullptr;
return false;
}

// convert x/y values back to distance/height values

Expand Down

0 comments on commit de9377a

Please sign in to comment.