Description
 

Sometimes, we may need to know how dense our point cloud is, that is, how many points we have per some area unit. With two key transofrmers and a few supporting transformers, it is an easy task with FME 2011.


First, we tile the source point cloud with Tiler - tile size defines area unit we will be using. The second key transformer, PointCloudPropertyExtractor, caclulates how many points fall into each tile.


After that, we force each tile to 3D, where the number of point is used as elevation, and finally, we rasterize 3D tiles to a numeric raster, which may look as, for example, the picture on the right side. This is a density raster of a point cloud (on the left) we thankfully received from the City of Gävle:


Point Cloud    


Density Raster



The roundish areas on the image appear due to the type of the scanner, which is constantly rotating.


In fact, the workflow is a little bit more complex. I used two Tilers, because if the number of small tiles exceeds a few hundreds of thousands, FME can run out of memory. For example, a 5000*5000 meters point cloud cannot be tiled to 1*1 m squares on our regular office machine.



Two-tier approach, where we first tile the point cloud into bigger tiles, and then tile each bigger tile to small, user-defined tiles, allows avoiding memory problems, or, at least push the limit to much higher numbers. Still, with more tiles, processing time grows exponentially. Check the following table:


Pixel Size    Translation time
    4               0'20"
    3               0'30"
    2               0'40"
    1               1'30"
    0.5             13'30"


With a few extra transformers, we can make the raster a bit more descriptive. RasterExpressionEvaluator with PaletteAdder allow to classify the UINT32 into an indexed raster, where each color signifies a certain range of points within a pixel:

This is how the expression in RasterExpressionEvaluator may look like:

if(A[0]==0,0,if(A[0]<5,1,if(A[0]<10,2,if(A[0]<25,3,if(A[0]<50,4,if(A[0]<100,5,if(A[0]<500,6,7)))))))

which means, that if the source pixel value is 0, keep it 0, if it is less than 5, make it 1, if it is less than 10, make it 2 and so on:



The output allows easily distinguish between different density zones on the raster:


Color Coded Raster   


Legend
(number of points per pixel)
   


Based on some internal tests, we may conclude that making such density rasters makes more sense for terrestrial LiDAR dataset.

The custom transformer is available for use in the attached file - PointCloudDensityCalculator.fmx