Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Fix elevation profile on virtual point cloud #59073

Merged
merged 2 commits into from
Oct 24, 2024
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
130 changes: 81 additions & 49 deletions src/core/pointcloud/qgspointcloudlayerprofilegenerator.cpp
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
Loading