Filter_LiDAR_DEM.sml

# Filter_LiDAR_DEM.sml - 	remove non-ground features from a LiDAR 
#													surface.
#
# This script filteres out non ground pixels (buildings, vegetation
# etc) from LiDAR surfaces (eg from LAS_to_Surface.sml)
#
#
# Method:
# -------
#
# 1.	A TIN is generated by finding the lowest pixels in a moving 
#			window, where the size of the moving window (horizontal) is 
#			set to be larger than the features to be removed. 
#
#	2.	A temporary surface raster is generated from the TIN.
#
#	3.	The slope of the TIN is calculated
#
#	4.	For each pixel where the slope does not exceed "max_slope"  
#			but where the DEM surface exceeds the TIN surface by more than
#			"vertical" make the pixel null. Where the DEM exceeds max_slope,
#			the surface is not filtered.
#
#	Best results are acheived by using the filter iteratively with
# successively smaller values for vertical and max_slope and each time
# reinstating pixels that are close to the gap filled output from the
# previous filter:
#
# Filter_LiDAR_DEM -> DEM_Interp -> Reinstate -> Filter_LiDAR_DEM etc.
#
#
# Input:  LiDAR surface (output from LAS_to_Raster.sml)
#	------
#
# Output: Filtered raster surface "non ground" features are NULL.
# -------
#
# Set the following paramaters to adjust the filtering.
#
numeric horizontal = 30; # Size of moving window for finding TIN nodes
numeric vertical = 2; 	# Max height above the DEM for ground features
numeric max_slope = 10;	# Maximum slope to filter.
#
# Terranean Mapping Technologies, Brisbane, Australia: 
# 20/12/2008
#####################################################################
clear();
class TIN T;
class RVC_RASTER DEM,Filtered,Temp, Slope;
numeric line,col, lines, cols, zmin, lstart, lfin, i, cstart, cfin, j, xmin, ymin;
GetInputRaster(DEM);
string rvc$ = GetObjectFileName(DEM);
string name$ = GetObjectName(rvc$,GetObjectNumber(DEM));
numeric nullValue = NullValue(DEM);
numeric cellSize = ColScale(DEM);
# Create Rasters
print("Creating temp rasters...");
string desc$ = "filter horizontal: "+NumToStr(max_slope)+", max_slope:"+NumToStr(vertical)+", vertical "+NumToStr(vertical);
CreateRaster(Filtered,rvc$,name$+"_flt",desc$,NumLins(DEM),NumCols(DEM),"32-bit float");
CreateRaster(Slope,rvc$,name$+"_slope",desc$,NumLins(DEM),NumCols(DEM),"32-bit float");
CreateTempRaster(Temp,NumLins(DEM),NumCols(DEM),"32-bit float");
CreateTIN(T,rvc$,"TIN","");
CopySubobjects(DEM,T,"georef");
CopySubobjects(DEM,Filtered,"georef");
CopySubobjects(DEM,Slope,"georef");
lines = NumLins(DEM);
cols = NumCols(DEM);
print("lines,cols",lines,cols);
print("pixel size",cellSize);
# Calulate spacings
numeric winsize = ceil(horizontal / (cellSize * 3)) * 3;
print("Rounded horizontal distance from",horizontal,"to",winsize*cellSize);
print("winsize",winsize);
numeric max_nodes = 2 * lines * cols / winsize;
array numeric xarr[max_nodes*2], yarr[max_nodes*2], zarr[max_nodes*2];
print("max_nodes",max_nodes);
numeric npts = 0;
SetStatusMessage("Finding nodes...");
# Find TIN nodes
for line = 1 to lines step winsize for col = 1 to cols step winsize
{	zmin = 9999;
	lstart = line - winsize/2;
	lfin = line + winsize/2;
	if (lstart < 1) lstart = 1;
	if (lfin > lines) lfin = lines;
	for  i = lstart to lfin
	{ 	cstart = col - winsize/2;
	 	cfin = col + winsize/2;
		if (cstart < 1) cstart = 1;
		if (cfin > cols) cfin = cols;
		for j = cstart to cfin
		{	if (DEM[i,j] < zmin) then
			{	zmin = DEM[i,j];
				xmin = j;
				ymin = i;
			}
		}
	}
	if (zmin <> 9999) then
	{	npts ++;
		xarr[npts] = xmin;
		yarr[npts] = ymin;
		zarr[npts] = zmin;
	}
}
# generate TIN
SetStatusMessage("Creating TIN,"+NumToStr(npts)+" nodes.");
TINCreateFromNodes(T,npts,xarr,yarr,zarr,1.0,0.0,0.0001);
SurfaceFitTINLinear(T,Temp,"","",DEM);
# Calculate Slope of TIN faces
foreach Temp Slope = FocalSlope(Temp) * deg;
# Search for non-ground features
SetStatusMessage("Filtering surface...");
foreach DEM if ( (DEM - Temp < vertical ) or (Slope > max_slope))   Filtered  = DEM;
SetStatusMessage("Finishing up...");
SetNull(Filtered,0);
CreateHistogram(Filtered);
CreatePyramid(Filtered);
SetStatusMessage("done.");
print("done.");
#eof