Skip to content

Commit 9a81fa7

Browse files
NathanW2NathanW2
committedMar 17, 2013
Merge pull request #461 from nyalldawson/heatmap_fixes
Massive speedup to heatmap plugin (fix #6756, fix #6691 and fix #6692)
·
ltr-3_44archive/master_2
2 parents 8966e9e + 239479a commit 9a81fa7

File tree

4 files changed

+86
-56
lines changed

4 files changed

+86
-56
lines changed
 

‎src/plugins/heatmap/heatmap.cpp

Lines changed: 66 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -107,8 +107,11 @@ void Heatmap::run()
107107
QgsRectangle myBBox = d.bbox();
108108
int columns = d.columns();
109109
int rows = d.rows();
110-
float cellsize = d.cellSizeX(); // or d.cellSizeY(); both have the same value
111-
float myDecay = d.decayRatio();
110+
double cellsize = d.cellSizeX(); // or d.cellSizeY(); both have the same value
111+
double myDecay = d.decayRatio();
112+
113+
// Start working on the input vector
114+
QgsVectorLayer* inputLayer = d.inputVectorLayer();
112115

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

130135
GDALRasterBand *poBand;
131136
poBand = emptyDataset->GetRasterBand( 1 );
@@ -153,23 +158,39 @@ void Heatmap::run()
153158
return;
154159
}
155160
poBand = heatmapDS->GetRasterBand( 1 );
156-
// Start working on the input vector
157-
QgsVectorLayer* inputLayer = d.inputVectorLayer();
158161

159162
QgsAttributeList myAttrList;
160163
int rField = 0;
161164
int wField = 0;
165+
166+
// Handle different radius options
167+
double radius;
168+
double radiusToMapUnits = 1;
169+
int myBuffer;
162170
if ( d.variableRadius() )
163171
{
164172
rField = d.radiusField();
165173
myAttrList.append( rField );
166174
QgsDebugMsg( QString( "Radius Field index received: %1" ).arg( rField ) );
175+
176+
// If not using map units, then calculate a conversion factor to convert the radii to map units
177+
if ( d.radiusUnit() == HeatmapGui::Meters )
178+
{
179+
radiusToMapUnits = mapUnitsOf( 1, inputLayer->crs() );
180+
}
181+
}
182+
else
183+
{
184+
radius = d.radius(); // radius returned by d.radius() is already in map units
185+
myBuffer = bufferSize( radius, cellsize );
167186
}
187+
168188
if ( d.weighted() )
169189
{
170190
wField = d.weightField();
171191
myAttrList.append( wField );
172192
}
193+
173194
// This might have attributes or mightnot have attibutes at all
174195
// based on the variableRadius() and weighted()
175196
QgsFeatureIterator fit = inputLayer->getFeatures( QgsFeatureRequest().setSubsetOfAttributes( myAttrList ) );
@@ -201,26 +222,14 @@ void Heatmap::run()
201222
{
202223
continue;
203224
}
204-
float radius;
225+
226+
// If radius is variable then fetch it and calculate new pixel buffer size
205227
if ( d.variableRadius() )
206228
{
207-
radius = myFeature.attribute( rField ).toFloat();
208-
}
209-
else
210-
{
211-
radius = d.radius();
212-
}
213-
//convert the radius to map units if it is in meters
214-
if ( d.radiusUnit() == HeatmapGui::Meters )
215-
{
216-
radius = mapUnitsOf( radius, inputLayer->crs() );
217-
}
218-
// convert radius in map units to pixel count
219-
int myBuffer = radius / cellsize;
220-
if ( radius - ( cellsize * myBuffer ) > 0.5 )
221-
{
222-
++myBuffer;
229+
radius = myFeature.attribute( rField ).toDouble() * radiusToMapUnits;
230+
myBuffer = bufferSize( radius, cellsize );
223231
}
232+
224233
int blockSize = 2 * myBuffer + 1; //Block SIDE would be more appropriate
225234
// calculate the pixel position
226235
unsigned int xPosition, yPosition;
@@ -232,18 +241,25 @@ void Heatmap::run()
232241
poBand->RasterIO( GF_Read, xPosition, yPosition, blockSize, blockSize,
233242
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );
234243

235-
float weight = 1.0;
244+
double weight = 1.0;
236245
if ( d.weighted() )
237246
{
238-
weight = myFeature.attribute( wField ).toFloat();
247+
weight = myFeature.attribute( wField ).toDouble();
239248
}
240249

241250
for ( int xp = 0; xp <= myBuffer; xp++ )
242251
{
243252
for ( int yp = 0; yp <= myBuffer; yp++ )
244253
{
245-
float distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
246-
float pixelValue = weight * ( 1 - (( 1 - myDecay ) * distance / myBuffer ) );
254+
double distance = sqrt( pow( xp, 2.0 ) + pow( yp, 2.0 ) );
255+
256+
// is pixel outside search bandwidth of feature?
257+
if ( distance > myBuffer )
258+
{
259+
continue;
260+
}
261+
262+
double pixelValue = weight * ( 1 - (( 1 - myDecay ) * distance / myBuffer ) );
247263

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

258-
if ( distance <= myBuffer )
274+
int pos[4];
275+
pos[0] = ( myBuffer + xp ) * blockSize + ( myBuffer + yp );
276+
pos[1] = ( myBuffer + xp ) * blockSize + ( myBuffer - yp );
277+
pos[2] = ( myBuffer - xp ) * blockSize + ( myBuffer + yp );
278+
pos[3] = ( myBuffer - xp ) * blockSize + ( myBuffer - yp );
279+
for ( int p = 0; p < 4; p++ )
259280
{
260-
int pos[4];
261-
pos[0] = ( myBuffer + xp ) * blockSize + ( myBuffer + yp );
262-
pos[1] = ( myBuffer + xp ) * blockSize + ( myBuffer - yp );
263-
pos[2] = ( myBuffer - xp ) * blockSize + ( myBuffer + yp );
264-
pos[3] = ( myBuffer - xp ) * blockSize + ( myBuffer - yp );
265-
for ( int p = 0; p < 4; p++ )
281+
if ( dataBuffer[ pos[p] ] == NO_DATA )
266282
{
267-
if ( dataBuffer[ pos[p] ] == NO_DATA )
268-
{
269-
dataBuffer[ pos[p] ] = 0;
270-
}
271-
dataBuffer[ pos[p] ] += pixelValue;
283+
dataBuffer[ pos[p] ] = 0;
272284
}
285+
dataBuffer[ pos[p] ] += pixelValue;
273286
}
274287
}
275288
}
@@ -278,7 +291,7 @@ void Heatmap::run()
278291
dataBuffer, blockSize, blockSize, GDT_Float32, 0, 0 );
279292
CPLFree( dataBuffer );
280293
}
281-
//Finally close the dataset
294+
// Finally close the dataset
282295
GDALClose(( GDALDatasetH ) heatmapDS );
283296

284297
// Open the file in QGIS window
@@ -291,7 +304,8 @@ void Heatmap::run()
291304
* Local functions
292305
*
293306
*/
294-
float Heatmap::mapUnitsOf( float meters, QgsCoordinateReferenceSystem layerCrs )
307+
308+
double Heatmap::mapUnitsOf( double meters, QgsCoordinateReferenceSystem layerCrs )
295309
{
296310
// Worker to transform metres input to mapunits
297311
QgsDistanceArea da;
@@ -304,6 +318,19 @@ float Heatmap::mapUnitsOf( float meters, QgsCoordinateReferenceSystem layerCrs )
304318
return meters / da.measureLine( QgsPoint( 0.0, 0.0 ), QgsPoint( 0.0, 1.0 ) );
305319
}
306320

321+
int Heatmap::bufferSize( double radius, double cellsize )
322+
{
323+
// Calculate the buffer size in pixels
324+
325+
int buffer = radius / cellsize;
326+
if ( radius - ( cellsize * buffer ) > 0.5 )
327+
{
328+
++buffer;
329+
}
330+
return buffer;
331+
}
332+
333+
307334
// Unload the plugin by cleaning up the GUI
308335
void Heatmap::unload()
309336
{

‎src/plugins/heatmap/heatmap.h

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -82,7 +82,10 @@ class Heatmap: public QObject, public QgisPlugin
8282

8383
private:
8484
//! Worker to convert meters to map units
85-
float mapUnitsOf( float meters, QgsCoordinateReferenceSystem crs );
85+
double mapUnitsOf( double meters, QgsCoordinateReferenceSystem layerCrs );
86+
//! Worker to calculate buffer size in pixels
87+
int bufferSize( double radius, double cellsize );
88+
8689
// MANDATORY PLUGIN PROPERTY DECLARATIONS .....
8790

8891
int mPluginType;

‎src/plugins/heatmap/heatmapgui.cpp

Lines changed: 10 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -170,7 +170,7 @@ void HeatmapGui::on_columnLineEdit_editingFinished()
170170

171171
void HeatmapGui::on_cellXLineEdit_editingFinished()
172172
{
173-
mXcellsize = cellXLineEdit->text().toFloat();
173+
mXcellsize = cellXLineEdit->text().toDouble();
174174
mYcellsize = mXcellsize;
175175
mRows = mBBox.height() / mYcellsize + 1;
176176
mColumns = mBBox.width() / mXcellsize + 1;
@@ -180,7 +180,7 @@ void HeatmapGui::on_cellXLineEdit_editingFinished()
180180

181181
void HeatmapGui::on_cellYLineEdit_editingFinished()
182182
{
183-
mYcellsize = cellYLineEdit->text().toFloat();
183+
mYcellsize = cellYLineEdit->text().toDouble();
184184
mXcellsize = mYcellsize;
185185
mRows = mBBox.height() / mYcellsize + 1;
186186
mColumns = mBBox.width() / mXcellsize + 1;
@@ -281,10 +281,10 @@ void HeatmapGui::updateBBox()
281281
mBBox = inputLayer->extent();
282282
QgsCoordinateReferenceSystem layerCrs = inputLayer->crs();
283283

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

289289
if ( radiusFieldUnitCombo->currentIndex() == HeatmapGui::Meters )
290290
{
@@ -297,7 +297,7 @@ void HeatmapGui::updateBBox()
297297
}
298298
else
299299
{
300-
float radiusValue = mBufferLineEdit->text().toFloat();
300+
double radiusValue = mBufferLineEdit->text().toDouble();
301301
if ( mRadiusUnitCombo->currentIndex() == HeatmapGui::Meters )
302302
{
303303
radiusInMapUnits = mapUnitsOf( radiusValue, layerCrs );
@@ -325,7 +325,7 @@ void HeatmapGui::updateBBox()
325325
}
326326
}
327327

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

359-
float HeatmapGui::radius()
359+
double HeatmapGui::radius()
360360
{
361-
float radius = mBufferLineEdit->text().toInt();
361+
double radius = mBufferLineEdit->text().toDouble();
362362
if ( mRadiusUnitCombo->currentIndex() == HeatmapGui::Meters )
363363
{
364364
radius = mapUnitsOf( radius, inputVectorLayer()->crs() );
@@ -375,9 +375,9 @@ int HeatmapGui::radiusUnit()
375375
return mRadiusUnitCombo->currentIndex();
376376
}
377377

378-
float HeatmapGui::decayRatio()
378+
double HeatmapGui::decayRatio()
379379
{
380-
return mDecayLineEdit->text().toFloat();
380+
return mDecayLineEdit->text().toDouble();
381381
}
382382

383383
int HeatmapGui::radiusField()

‎src/plugins/heatmap/heatmapgui.h

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -43,13 +43,13 @@ class HeatmapGui : public QDialog, private Ui::HeatmapGuiBase
4343
bool variableRadius();
4444

4545
/** Returns the fixed radius value */
46-
float radius();
46+
double radius();
4747

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

5151
/** Return the decay ratio */
52-
float decayRatio();
52+
double decayRatio();
5353

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

7575
/** Returns the cell size X value */
76-
float cellSizeX() { return mXcellsize; }
76+
double cellSizeX() { return mXcellsize; }
7777

7878
/** Returns the cell size Y valuue */
79-
float cellSizeY() { return mYcellsize; }
79+
double cellSizeY() { return mYcellsize; }
8080

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

8787
// bbox of layer for lineedit changes
8888
QgsRectangle mBBox;
89-
float mXcellsize, mYcellsize;
89+
double mXcellsize, mYcellsize;
9090
int mRows, mColumns;
9191

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

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

107107
private slots:
108108
void on_mButtonBox_accepted();

0 commit comments

Comments
 (0)
Please sign in to comment.