Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
Merge pull request #461 from nyalldawson/heatmap_fixes
Massive speedup to heatmap plugin (fix #6756, fix #6691 and fix #6692)
  • Loading branch information
NathanW2 committed Mar 17, 2013
2 parents 8966e9e + 239479a commit 9a81fa7
Show file tree
Hide file tree
Showing 4 changed files with 86 additions and 56 deletions.
105 changes: 66 additions & 39 deletions src/plugins/heatmap/heatmap.cpp
Expand Up @@ -107,8 +107,11 @@ void Heatmap::run()
QgsRectangle myBBox = d.bbox();
int columns = d.columns();
int rows = d.rows();
float cellsize = d.cellSizeX(); // or d.cellSizeY(); both have the same value
float myDecay = d.decayRatio();
double cellsize = d.cellSizeX(); // or d.cellSizeY(); both have the same value
double myDecay = d.decayRatio();

// Start working on the input vector
QgsVectorLayer* inputLayer = d.inputVectorLayer();

// Getting the rasterdataset in place
GDALAllRegister();
Expand All @@ -126,6 +129,8 @@ void Heatmap::run()
double geoTransform[6] = { myBBox.xMinimum(), cellsize, 0, myBBox.yMinimum(), 0, cellsize };
emptyDataset = myDriver->Create( d.outputFilename().toUtf8(), columns, rows, 1, GDT_Float32, NULL );
emptyDataset->SetGeoTransform( geoTransform );
// Set the projection on the raster destination to match the input layer
emptyDataset->SetProjection( inputLayer->crs().toWkt().toLocal8Bit().data() );

GDALRasterBand *poBand;
poBand = emptyDataset->GetRasterBand( 1 );
Expand Down Expand Up @@ -153,23 +158,39 @@ void Heatmap::run()
return;
}
poBand = heatmapDS->GetRasterBand( 1 );
// Start working on the input vector
QgsVectorLayer* inputLayer = d.inputVectorLayer();

QgsAttributeList myAttrList;
int rField = 0;
int wField = 0;

// Handle different radius options
double radius;
double radiusToMapUnits = 1;
int myBuffer;
if ( d.variableRadius() )
{
rField = d.radiusField();
myAttrList.append( rField );
QgsDebugMsg( QString( "Radius Field index received: %1" ).arg( rField ) );

// If not using map units, then calculate a conversion factor to convert the radii to map units
if ( d.radiusUnit() == HeatmapGui::Meters )
{
radiusToMapUnits = mapUnitsOf( 1, inputLayer->crs() );
}
}
else
{
radius = d.radius(); // radius returned by d.radius() is already in map units
myBuffer = bufferSize( radius, cellsize );
}

if ( d.weighted() )
{
wField = d.weightField();
myAttrList.append( wField );
}

// This might have attributes or mightnot have attibutes at all
// based on the variableRadius() and weighted()
QgsFeatureIterator fit = inputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( myAttrList ) );
Expand Down Expand Up @@ -201,26 +222,14 @@ void Heatmap::run()
{
continue;
}
float radius;

// If radius is variable then fetch it and calculate new pixel buffer size
if ( d.variableRadius() )
{
radius = myFeature.attribute( rField ).toFloat();
}
else
{
radius = d.radius();
}
//convert the radius to map units if it is in meters
if ( d.radiusUnit() == HeatmapGui::Meters )
{
radius = mapUnitsOf( radius, inputLayer->crs() );
}
// convert radius in map units to pixel count
int myBuffer = radius / cellsize;
if ( radius - ( cellsize * myBuffer ) > 0.5 )
{
++myBuffer;
radius = myFeature.attribute( rField ).toDouble() * radiusToMapUnits;
myBuffer = bufferSize( radius, cellsize );
}

int blockSize = 2 * myBuffer + 1; //Block SIDE would be more appropriate
// calculate the pixel position
unsigned int xPosition, yPosition;
Expand All @@ -232,18 +241,25 @@ void Heatmap::run()
poBand->RasterIO( GF_Read, xPosition, yPosition, blockSize, blockSize,
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );

float weight = 1.0;
double weight = 1.0;
if ( d.weighted() )
{
weight = myFeature.attribute( wField ).toFloat();
weight = myFeature.attribute( wField ).toDouble();
}

for ( int xp = 0; xp <= myBuffer; xp++ )
{
for ( int yp = 0; yp <= myBuffer; yp++ )
{
float distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
float pixelValue = weight * ( 1 - (( 1 - myDecay ) * distance / myBuffer ) );
double distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );

// is pixel outside search bandwidth of feature?
if ( distance > myBuffer )
{
continue;
}

double pixelValue = weight * ( 1 - (( 1 - myDecay ) * distance / myBuffer ) );

// clearing anamolies along the axes
if ( xp == 0 && yp == 0 )
Expand All @@ -255,21 +271,18 @@ void Heatmap::run()
pixelValue /= 2;
}

if ( distance <= myBuffer )
int pos[4];
pos[0] = ( myBuffer + xp ) * blockSize + ( myBuffer + yp );
pos[1] = ( myBuffer + xp ) * blockSize + ( myBuffer - yp );
pos[2] = ( myBuffer - xp ) * blockSize + ( myBuffer + yp );
pos[3] = ( myBuffer - xp ) * blockSize + ( myBuffer - yp );
for ( int p = 0; p < 4; p++ )
{
int pos[4];
pos[0] = ( myBuffer + xp ) * blockSize + ( myBuffer + yp );
pos[1] = ( myBuffer + xp ) * blockSize + ( myBuffer - yp );
pos[2] = ( myBuffer - xp ) * blockSize + ( myBuffer + yp );
pos[3] = ( myBuffer - xp ) * blockSize + ( myBuffer - yp );
for ( int p = 0; p < 4; p++ )
if ( dataBuffer[ pos[p] ] == NO_DATA )
{
if ( dataBuffer[ pos[p] ] == NO_DATA )
{
dataBuffer[ pos[p] ] = 0;
}
dataBuffer[ pos[p] ] += pixelValue;
dataBuffer[ pos[p] ] = 0;
}
dataBuffer[ pos[p] ] += pixelValue;
}
}
}
Expand All @@ -278,7 +291,7 @@ void Heatmap::run()
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );
CPLFree( dataBuffer );
}
//Finally close the dataset
// Finally close the dataset
GDALClose(( GDALDatasetH ) heatmapDS );

// Open the file in QGIS window
Expand All @@ -291,7 +304,8 @@ void Heatmap::run()
* Local functions
*
*/
float Heatmap::mapUnitsOf( float meters, QgsCoordinateReferenceSystem layerCrs )

double Heatmap::mapUnitsOf( double meters, QgsCoordinateReferenceSystem layerCrs )
{
// Worker to transform metres input to mapunits
QgsDistanceArea da;
Expand All @@ -304,6 +318,19 @@ float Heatmap::mapUnitsOf( float meters, QgsCoordinateReferenceSystem layerCrs )
return meters / da.measureLine( QgsPoint( 0.0, 0.0 ), QgsPoint( 0.0, 1.0 ) );
}

int Heatmap::bufferSize( double radius, double cellsize )
{
// Calculate the buffer size in pixels

int buffer = radius / cellsize;
if ( radius - ( cellsize * buffer ) > 0.5 )
{
++buffer;
}
return buffer;
}


// Unload the plugin by cleaning up the GUI
void Heatmap::unload()
{
Expand Down
5 changes: 4 additions & 1 deletion src/plugins/heatmap/heatmap.h
Expand Up @@ -82,7 +82,10 @@ class Heatmap: public QObject, public QgisPlugin

private:
//! Worker to convert meters to map units
float mapUnitsOf( float meters, QgsCoordinateReferenceSystem crs );
double mapUnitsOf( double meters, QgsCoordinateReferenceSystem layerCrs );
//! Worker to calculate buffer size in pixels
int bufferSize( double radius, double cellsize );

// MANDATORY PLUGIN PROPERTY DECLARATIONS .....

int mPluginType;
Expand Down
20 changes: 10 additions & 10 deletions src/plugins/heatmap/heatmapgui.cpp
Expand Up @@ -170,7 +170,7 @@ void HeatmapGui::on_columnLineEdit_editingFinished()

void HeatmapGui::on_cellXLineEdit_editingFinished()
{
mXcellsize = cellXLineEdit->text().toFloat();
mXcellsize = cellXLineEdit->text().toDouble();
mYcellsize = mXcellsize;
mRows = mBBox.height() / mYcellsize + 1;
mColumns = mBBox.width() / mXcellsize + 1;
Expand All @@ -180,7 +180,7 @@ void HeatmapGui::on_cellXLineEdit_editingFinished()

void HeatmapGui::on_cellYLineEdit_editingFinished()
{
mYcellsize = cellYLineEdit->text().toFloat();
mYcellsize = cellYLineEdit->text().toDouble();
mXcellsize = mYcellsize;
mRows = mBBox.height() / mYcellsize + 1;
mColumns = mBBox.width() / mXcellsize + 1;
Expand Down Expand Up @@ -281,10 +281,10 @@ void HeatmapGui::updateBBox()
mBBox = inputLayer->extent();
QgsCoordinateReferenceSystem layerCrs = inputLayer->crs();

float radiusInMapUnits = 0.0;
double radiusInMapUnits = 0.0;
if ( useRadius->isChecked() )
{
float maxInField = inputLayer->maximumValue( radiusFieldCombo->itemData( radiusFieldCombo->currentIndex() ).toInt() ).toFloat();
double maxInField = inputLayer->maximumValue( radiusFieldCombo->itemData( radiusFieldCombo->currentIndex() ).toInt() ).toDouble();

if ( radiusFieldUnitCombo->currentIndex() == HeatmapGui::Meters )
{
Expand All @@ -297,7 +297,7 @@ void HeatmapGui::updateBBox()
}
else
{
float radiusValue = mBufferLineEdit->text().toFloat();
double radiusValue = mBufferLineEdit->text().toDouble();
if ( mRadiusUnitCombo->currentIndex() == HeatmapGui::Meters )
{
radiusInMapUnits = mapUnitsOf( radiusValue, layerCrs );
Expand Down Expand Up @@ -325,7 +325,7 @@ void HeatmapGui::updateBBox()
}
}

float HeatmapGui::mapUnitsOf( float meters, QgsCoordinateReferenceSystem layerCrs )
double HeatmapGui::mapUnitsOf( double meters, QgsCoordinateReferenceSystem layerCrs )
{
// converter function to transform metres input to mapunits
// so that bounding box can be updated
Expand Down Expand Up @@ -356,9 +356,9 @@ bool HeatmapGui::variableRadius()
return useRadius->isChecked();
}

float HeatmapGui::radius()
double HeatmapGui::radius()
{
float radius = mBufferLineEdit->text().toInt();
double radius = mBufferLineEdit->text().toDouble();
if ( mRadiusUnitCombo->currentIndex() == HeatmapGui::Meters )
{
radius = mapUnitsOf( radius, inputVectorLayer()->crs() );
Expand All @@ -375,9 +375,9 @@ int HeatmapGui::radiusUnit()
return mRadiusUnitCombo->currentIndex();
}

float HeatmapGui::decayRatio()
double HeatmapGui::decayRatio()
{
return mDecayLineEdit->text().toFloat();
return mDecayLineEdit->text().toDouble();
}

int HeatmapGui::radiusField()
Expand Down
12 changes: 6 additions & 6 deletions src/plugins/heatmap/heatmapgui.h
Expand Up @@ -43,13 +43,13 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase
bool variableRadius();

/** Returns the fixed radius value */
float radius();
double radius();

/** Return the radius Unit (meters/map units) */
int radiusUnit();

/** Return the decay ratio */
float decayRatio();
double decayRatio();

/** Return the attribute field for variable radius */
int radiusField();
Expand All @@ -73,10 +73,10 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase
int columns() { return mColumns; }

/** Returns the cell size X value */
float cellSizeX() { return mXcellsize; }
double cellSizeX() { return mXcellsize; }

/** Returns the cell size Y valuue */
float cellSizeY() { return mYcellsize; }
double cellSizeY() { return mYcellsize; }

/** Return the BBox */
QgsRectangle bbox() { return mBBox; }
Expand All @@ -86,7 +86,7 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase

// bbox of layer for lineedit changes
QgsRectangle mBBox;
float mXcellsize, mYcellsize;
double mXcellsize, mYcellsize;
int mRows, mColumns;

/** Function to check wether all constrains are satisfied and enable the OK button */
Expand All @@ -102,7 +102,7 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase
void updateSize();

/** Convert Maters value to the corresponding map units based on Layer projection */
float mapUnitsOf( float meters, QgsCoordinateReferenceSystem layerCrs );
double mapUnitsOf( double meters, QgsCoordinateReferenceSystem layerCrs );

private slots:
void on_mButtonBox_accepted();
Expand Down

0 comments on commit 9a81fa7

Please sign in to comment.