interp.sml

# INTERP.SML - Interpolate gaps in a surface raster.
#
# This script is designed to fill holes in a DEM generated from high
# density (relative to cell size) LiDAR and gaps left after raster 
#  filtering non-ground features from a DEM.
#
# The algorithm is a two stage process:
#
#	1st, small gaps are filled by averaging, 
#	2nd, large gaps are filled by surface fitting.
#
# Outline:
# --------
# 1. Copy input raster (DEM) to output raster (Interp), filling small
#    gaps with a 3x3 focal mean.
# 2. Cells adjacent to unfilled gaps are loaded into arrays. (xarr, 
#    yarr & zarr) 
# 3. A temporary TIN is generated from the gap edge coordinates. (T)
# 4. A temporary Raster surface is fitted to the TIN. (Temp1)
# 5. The Raster surface is resampled to match the input DEM. (Temp2)
# 6. Large gaps are filled by transferring cell values from TIN 
#    raster Temp2 to gaps in Interp.
#
# Input:  DEM with gaps
# Output: DEM without gaps - output raster is created in same rvc 
#	  file "_interp" is appended to the object name 
#
# Constraints:  The required size of arrays xarr, yarr, & zarr can't be 
#		known in advance - possible overflow error.
#
#too many nodes, increase "max_nodes" on the next line:
numeric max_nodes = 1000000;
# Gaps in input surface must be set NULL.
# Any valid surface values of 0.00 will be reinterpolated.
#
# Terranean Mapping Technologies, Brisbane, Australia
# 20/11/2008
#####################################################################
clear();
class TIN T;
class RASTER DEM,Interp,Temp1, Temp2;
array numeric xarr[max_nodes], yarr[max_nodes], zarr[max_nodes];
numeric line,col;
#Get Input raster
GetInputRaster(DEM);
string rvc$ = GetObjectFileName(DEM);
string name$ = GetObjectName(rvc$,GetObjectNumber(DEM));
numeric nullValue = NullValue(DEM);
#Create Temp Rasters
print("Creating temp rasters...");
string desc$ = "Interp "+DEM.$Info.Desc;
CreateRaster(Interp,rvc$,name$+"_interp",desc$,NumLins(DEM),NumCols(DEM),"32-bit float");
CreateTempRaster(Temp1,NumLins(DEM),NumCols(DEM),"32-bit float");
CreateTempRaster(Temp2,NumLins(DEM),NumCols(DEM),"32-bit float");
CreateTIN(T,rvc$,"TIN","");
CopySubobjects(DEM,T,"georef");
CopySubobjects(DEM,Interp,"georef");
#Fill small holes using Focal mean 3x3 window.
print("Filling small holes...");
foreach DEM
	if (DEM <> nullValue) then
		Interp = DEM;
	else
		Interp = FocalMean(DEM,3,3) ;
#Fill big holes using TIN method
print("Finding big holes...");
numeric npts = 0;
foreach Interp[line,col] if (Interp <> 0) then
	if (FocalMin(Interp,3,3) == 0) then
	{	npts ++;
		if (npts > max_nodes) then
		{	string msg$ = "\n\nArray size 'max_nodes' ("+NumToStr(max_nodes)+") exceeded\n\n Increase max_nodes and try again.\n\n";
			PopupMessage(msg$);
			Exit();
		}
		xarr[npts] = col;
		yarr[npts] = line;
		zarr[npts] = Interp[line,col];
	}
if (npts == 0) then
{	PopupMessage("\n\n    Didn't find any holes. Check null value.   \n\n" );
	Exit();
}
print("Creating TIN",npts,"nodes...");
TINCreateFromNodes(T,npts,xarr,yarr,zarr,1.0,0.0,0.0001);
SurfaceFitTINLinear(T,Temp1);
print("Resampling TIN raster...");
ResampleRasterToMatch(Temp1,Interp,Temp2,"affine","bilinear");
DeleteObject(rvc$,GetObjectNumber(T));
print("Filling big holes...");
foreach Interp[line,col]  if (Interp == 0) Interp = Temp2;
#Finish up
print("Finishing up...");
SetNull(Interp,0);
CreateHistogram(Interp);
CreatePyramid(Interp);
print("done.");
#eof