### surfcurve.sml

### Computes curvatures (profile = vertical and plan = horizontal) 
### for each DEM raster cell.  Profile curvature is positive for
### convex upward surfaces; plan curvature is positive for convex
### downhill surfaces. 

### Curvatures are computed for each DEM cell using a 3 x 3 kernel
### by computing a local best-fit quadratic surface for the nine-cell 
### neighborhood.  The equation of this surface can be written:

### z = a*x^2 + b*y^2 + c*xy + d*x + ee*y + f

### The best-fit surface is not constrained to exactly fit any of the
### DEM values in the neighborhood; it provides a smoothed surface
### that minimizes the impact of elevation errors in the DEM. 

### Algorithm modified from:
### Evans, I.S. (1979) An integrated system of terrain analysis and slope mapping. 
### Final report on grant DA-ERO-591-73-G0040, University of Durham, England.
### http://www.soi.city.ac.uk/~jwo/phd/04param.php3

### See also Schmidt, J., Evans, I.S., and Brinkmann, J. (2003), Comparison of
### polynomial models for land surface curvature calculation.  International
### Journal of Geographical Information Science, v.17, no.8, p. 797-814.

#### Global variables #########################################
 
numeric scale;
### With scale = 1, units of curvature are radians / meter.
### With scale = 100, units of curvature are radians / 100 m.
 
numeric a, b, c, d, ee, f;		# coefficients of quadratic equation
raster DEM;							# input elevation surface raster
raster Profile, Plan;			# output curvature rasters
numeric lin, col;					# line and column counters
numeric cellX, cellY;			# DEM cell sizes in column and line directions
numeric prfdenom, plndenom;	# values of denominators of the two curvature expressions

numeric doSmooth;					# flag for whether or not to used averaged DEM value
array numeric kernel[3,3];

# Procedure to read 3x3 kernel of DEM values
# surrounding current cell.  Notation of cells in kernel is as follows:

#	 ----------------
#	 |	z1 | z2 | z3 |
#	 ----------------
#	 |	z4	| z5 | z6 |
#   ----------------
#	 |	z7	| z8 | z9 |
#	 ----------------

proc readkernel (class raster K, numeric line, numeric colm){
	numeric z1 = K[line-1, colm-1];
	numeric z2 = K[line-1, colm];
	numeric z3 = K[line-1, colm+1];
	numeric z4 = K[line, colm-1];
	numeric z5 = K[line, colm];
	numeric z6 = K[line, colm+1];
	numeric z7 = K[line+1, colm-1];
	numeric z8 = K[line+1, colm];
	numeric z9 = K[line+1, colm+1];
	}  

# Function to check surrounding area of Kernel.  If no difference, 
# there is no need to do some computations.
# Also used as a safety against divide by zero.
func levelcheck() {
	if (z5 != z1) return true;
	if (z5 != z2) return true;
 	if (z5 != z3) return true;
	if (z5 != z4) return true;
	if (z5 != z6) return true;
	if (z5 != z7) return true;
	if (z5 != z8) return true; 
	if (z5 != z9) return true;
	return false;
}

# Function to check kernel values for null cells.
# Returns false if any cell is null. 
func nullcheck() {
	if (IsNull(z1) ) return false;
	if (IsNull(z2) ) return false;
	if (IsNull(z3) ) return false;
	if (IsNull(z4) ) return false;
	if (IsNull(z5) ) return false;
	if (IsNull(z6) ) return false;
	if (IsNull(z7) ) return false;
	if (IsNull(z8) ) return false;
	if (IsNull(z9) ) return false;
	return true;
	}

## procedure to compute coefficients for local quadratic surface
## dx and dy are line and column cell sizes
proc coeff(numeric dx, numeric dy) {
	a = (z1 + z3 + z4 + z6 + z7 + z9)/(6 * dx^2) - (z2 + z5 + z8)/(3 * dx^2);
	b = (z1 + z2 + z3 + z7 + z8 + z9)/(6 * dy^2) - (z4 + z5 + z6)/(3 * dy^2);
	c = (z3 + z7 - z1 - z9)/(4 * dx * dy);
	d = (z3 + z6 + z9 - z1 - z4 - z7)/(6 * dx);
	ee = (z1 + z2 + z3 - z7 - z8 -z9)/(6 * dy);
	}


proc makeCurv(class raster D) {

	for each D[lin,col] {

		if (lin != 1 && lin != D.$Info.NumLins && col != 1 && col != D.$Info.NumCols) {  ## don't process edge cells

			readkernel(D, lin, col);		# read 3x3 kernel values and check for null cells
			if ( nullcheck() ) {

				if( levelcheck() ) {		## don't process cell if 3x3 kernel is level

					coeff(cellX, cellY);		## call procedure to compute coefficients for quadratic surface

					# compute denominators for profile and plan curvature expressions to check for division by zero
					prfdenom = round(10^7 * ( (sqr(ee) + sqr(d) ) * ( (1 + sqr(d) + sqr(ee) ) ^ 1.5)) ) * 10^-7;
					plndenom = round(10^7 * ( (sqr(ee) + sqr(d) )^1.5) ) * 10^-7;

					if (prfdenom == 0 || plndenom == 0) {		## avoid division by zero (can occur for summits, depressions, and saddles)

						if ( (a > 0 && b > 0) || (a < 0 && b < 0) ) { 
							Profile[lin, col] = - (a + b) * scale;				## case for summit or depression
							Plan[lin, col] = -32768;
							}
						else {
							Profile[lin,col] = 0;					## case for saddles
							Plan[lin, col] = 0;
							}
						}		# end if (XXXdenom == 0)
					else {		# curvature computation for typical cells
						Profile[lin, col] = scale * -2 * ( a * sqr(d) + b * sqr(ee) + c * d * ee) / prfdenom;
						Plan[lin, col] = scale * -2 * (b * sqr(d) + a * sqr(ee) - c * d * ee) / plndenom;
						}
					} 
				else {
					Plan[lin, col] = -32768;			## case for flat kernel area
					Profile[lin, col] = 0;
					}
				} # end if ( levelcheck() )

			else {
				# case for kernels including null cells, not enough data to compute curvature
				Plan[lin,col] =-32768;
				Profile[lin,col] =-32768;
				}
			} # end if ( nullcheck() )
	else {
		# case for edge cells, not enough data to compute curvature
		Plan[lin,col] = -32768;
		Profile[lin,col] = -32768;
		}
	}
}


#############################################################
### Main program

GetInputRaster(DEM);
cellX = ColScale(DEM);
cellY = LinScale(DEM);

GetOutputRaster(Profile, DEM.$Info.NumLins, DEM.$Info.NumCols, "32-bit float");
GetOutputRaster(Plan, DEM.$Info.NumLins, DEM.$Info.NumCols, "32-bit float");
SetNull(Profile,-32768);
SetNull(Plan,-32768);

scale = PopupNum("Enter 1 to compute curvature in radians/m or \nenter 100 for curvature in radians/100 m", 1, 100, 0);

doSmooth = PopupYesNo("Use averaged DEM values for curvature?", 1);

if (doSmooth == 1) {
	class RASTER Temp;
	# pop up dialog window to get size of filter window for averaging DEM values
	local numeric filtsize = PopupNum("Size of averaging window in cells? \n(3, 5, 7, or 9)", 5, 3, 9, 0);

	if ( filtsize % 2 == 0 ) then		## check for even number entered as window size
		filtsize = filtsize + 1;

	CreateTempRaster(Temp, DEM.$Info.NumLins, DEM.$Info.NumCols, DEM.$Info.Type);

	## compute averaged DEM value for non-null cell, else assign null value to Temp raster
	numeric i, j;
	for i = 1 to DEM.$Info.NumLins  {
		for j = 1 to DEM.$Info.NumCols {
			if ( IsNull(DEM[i,j]) ) then
				Temp[i,j] = -32768;
			else
				Temp[i,j] = FocalMean(DEM, filtsize, filtsize, i, j);
			}
		}

	SetNull(Temp, -32768);

	makeCurv(Temp);
	CloseRaster(Temp);
	}	# end of if (doSmooth() )

else
	makeCurv(DEM);

for each Profile {
	if (Profile < -32768) then
		Profile = -32768;

	if (Plan < -32768) then
		Plan = -32768;
	}


SetNull(Profile, -32768);
SetNull(Plan, -32768);
CreateHistogram(Profile);
CreateHistogram(Plan);
CreatePyramid(Profile);
CreatePyramid(Plan);
CopySubobjects(DEM, Profile, "georef");
CopySubobjects(DEM, Plan, "georef");