######################################################################################################
#
# MorphContour.sml  (Stand Alone)
#
######################################################################################################
#
# Author:
#    Dave Breitwisch, MicroImages, Inc.
#    with contributions from Randy Smith, MicroImages, Inc.
#
# Started:
#    17 November 2005
#
# Updated:
#    18 July 2006
#
######################################################################################################
#
# Documentation:  MorphContour.sml
#
# Purpose:
#
#  The purpose of this Stand Alone SML script is to demonstrate the results of a surface-fitting
#  process from contour lines using morphological functions as described in the paper "An Image Space 
#  Algorithm for Morphological Contour Interpolation" by Barrett, Mortensen and Taylor.
#
#  The basic concept of this algorithm is to grow two known contour lines until their medial line is 
#  reached.  This line is then given the average value of the two lines and is added to the resulting 
#  image.  This new line is then used in the next iteration to find additional medial lines.  The 
#  process continues until there are no new cell values being added to the final image.
#
#  Some limitations of this algorithm include: 1) no filling of areas that are in between contour lines of 
#  the same value, such as ridges, saddles, peaks, or depressions; 2) there are also limitations based 
#  on the element structure used in the morphological functions, which lead to some unfilled cells if
#  the output DEM is set to be an integer data type.
#
######################################################################################################

# Function Prototypes
func MyErosion(class RASTER Rin, class RASTER Rtemp, numeric maskValue);
proc MyDilation(class RASTER Rin, class RASTER Rout, class RASTER Rmask, numeric maskValue, numeric minValue);
proc MyRasterXOR(class RASTER Rin1, class RASTER Rin2, class RASTER Rout);
proc MyMedialLine(class RASTER Rin1, class RASTER Rin2, class RASTER Rmask);
func MyAddToAccumulator(class RASTER Rin, class RASTER Rmask, class RASTER Rout, numeric intContValue);
proc Main();

# Global Variables
class VECTOR Vect;		# the input vector containing the 3D contours
class RASTER A;			# raster with initial contour z values transfered from input vector
class RASTER T, U, V, W, Temp;	# intermediate working rasters
class RASTER DEM;		# the final output DEM raster
numeric numlins, numcols;
array numeric element[3,3];
string table$, field$;		# names of line table and field containing elevation values

# Used as global rather then passing by reference since they are used in recursive functions
class RASTER FILL; #intermediate raster for filling holes
numeric value;
numeric currsize = 64;
numeric currpos = 0;
array numeric rowlist[currsize];
array numeric collist[currsize];

# Start process here
######################################################################################################
clear();
Main();

# Additional functions called
######################################################################################################


######################################################################################################
# This is the main function for the script.
proc Main() {

	numeric done = 1, numRuns = 0;
	# Create the Morphological element
	element[1,1] = null;
	element[1,2] = 0;
	element[1,3] = null;
	element[2,1] = 0;
	element[2,2] = 0;
	element[2,3] = 0;
	element[3,1] = null;
	element[3,2] = 0;
	element[3,3] = null;

	local class GEOREF vGeoref, rGeoref;
	local numeric dbErr;
	local class DATABASE vlineDB;		# vector line database
	local class DBTABLEINFO tInfo;	# class instance for selected line table
	local class DBFIELDINFO fInfo;	# class instance for selected field with elevation values

	local string outFilename$, outRastname$, outRastdesc$, tempFilename$;
	local string prompt$;
	local numeric rType;			# numeric flag for output data type for DEM
	local string rastType$;		# raster data type string
	local numeric cellSize;
	local numeric k, elev;

	# select the input vector object
	PopupMessage("Select a vector object containing the vector contours with elevations:");
	GetInputVector(Vect);
	vGeoref = GetLastUsedGeorefObject(Vect);
	if (vGeoref == 0) {
		CloseVector(Vect);
		PopupMessage("The vector you selected does not contain a valid georeference.  Exiting the script....");
		return;
	}

	# select the table and field in the line database that contains the elevation values
	vlineDB = OpenVectorLineDatabase(Vect);
	PopupMessage("Select line database table and field containing the contour line elevations:");
	dbErr = PopupSelectTableField(vlineDB, table$, field$);

	if (dbErr < 1) {		# exit script if user cancels the table selection dialog
		printf("No database table selected; exiting script.");
		Exit();
		}

	tInfo = DatabaseGetTableInfo(vlineDB, table$);
	fInfo = FieldGetInfoByName(tInfo, field$);

	# prompt for the output project file and name and description for the output DEM
	outFilename$ = GetOutputFileName("","Select or create a project file for the output DEM:", "rvc");
	outRastname$ = PopupString("Enter a name for the output DEM raster:", "DEM");
	outRastdesc$ = PopupString("Enter a description for the output DEM raster:", "");

	# prompt for and set the data type for the output DEM raster
	prompt$ = "Choose bit-depth of output elevation raster:\n" +
					"Press Yes to accept the default: 16-bit signed integer.\n" +
					"Press No to change to 32-bit floating point values.";
	rType = PopupYesNo(prompt$, 1);

	if (rType == 1) then
		rastType$ = "16-bit signed";
	else
		rastType$ = "32-bit float";

	# prompt for the cell size for the output DEM
	cellSize = PopupNum("Enter desired cell size of output raster in meters", 30, 0);


	# create temporary file for the contour raster R 
	tempFilename$ = CreateTempFileName();

	# create a blank contour raster A using the extents of the input vector and the specified cell size
	CreateRasterFromObject(Vect, A, tempFilename$, "RastCon", "", cellSize, cellSize, rastType$, vGeoref);
	rGeoref = GetLastUsedGeorefObject(A);
	numlins = NumLins(A);
	numcols = NumCols(A);

	# tranfer contour line elevation values to the contour raster
	for k = 1 to Vect.$Info.NumLines {
		elev = TableReadFieldNum(tInfo, field$, k);
		VectorElementToRaster(Vect, A, "line", k, elev, vGeoref, rGeoref);
		}

	# Use Max and Min values outside the boundaries of the contour line values -or- use the Max and Min values of the raster
	numeric max = GlobalMax(A) + 1;
	#numeric Max = Vect.$Info.MaxZ + 1;
	numeric min = GlobalMin(A);
	#numeric Min = Vect.$Info.MinZ - 1;

	# Create the necessary working objects
	CreateRaster(T, tempFilename$, "T", "T", numlins, numcols, rastType$);
	CreateRaster(U, tempFilename$, "U", "U", numlins, numcols, rastType$);
	CreateRaster(V, tempFilename$, "V", "V", numlins, numcols, "binary");
	CopyGeorefToObject(V, rGeoref);
	CreateRaster(W, tempFilename$, "W", "W", numlins, numcols, rastType$);
	CreateTempRaster(Temp, numlins, numcols, rastType$, 0); 

	while (done) {

		done = 0;
		numRuns++;

		# Filling all inter-contour cells in mask rasters with the maximum value 
		RasterCopy(A, W);
		RasterCopy(A, T);
		numeric i, j;
		for i = 1 to numlins {
			for j = 1 to numcols {
				if (A[i,j] == min) {
					W[i,j] = max;
					T[i,j] = max;
				}
			}
		}

		# use morphological functions to compute a set of medial lines for all contour pairs
		while ( MyErosion(W, Temp, max) > 0) { }
		MyDilation(W, U, T, max, min);
		MyRasterXOR(W, U, V);
		MyMedialLine(W, U, V);
		done = MyAddToAccumulator(W, V, A, min);

		print("Accumulator Count:", done, "\nRun Completed:", numRuns, "\n");

	}

	# close intermediate work rasters used in the morphological interpolation
	CloseRaster(T);
	CloseRaster(W);
	CloseRaster(U);
	CloseRaster(V);

	# Fill blank areas left after morphological interpolation
	CreateRaster(FILL, tempFilename$, "FILL", "FILL", numlins, numcols, rastType$);
	RasterCopy(A, FILL);
	# Set min values to maximum before filling by erosion
	for i = 1 to numlins {
		for j = 1 to numcols {
			if (FILL[i,j] == min) FILL[i,j] = max;
		}
	}
	while ( MyErosion(FILL, Temp, max) > 0 ) { } 

	CloseRaster(Temp);

	CreateRaster(DEM, outFilename$, outRastname$, outRastdesc$, numlins, numcols, rastType$);

	RasterCopy(FILL, DEM);
	CreatePyramid(DEM);
	CreateHistogram(DEM);
	CloseRaster(DEM);

	CloseRaster(A);
	CloseRaster(FILL);
}

######################################################################################################
# This function copies the input raster to the output raster (Rtemp) and morphologically
# erodes the output raster cells that have a starting value equal to the maskValue (inter-contour
# cells).  The function returns the number of cells that were eroded.
func MyErosion(class RASTER Rin, class RASTER Rtemp, numeric maskValue) {

	RasterCopy(Rin, Rtemp);

	local numeric r, c, er, ec, tempValue, newValue, count=0;
	for r = 1 to NumLins(Rin) {
		for c = 1 to NumCols(Rin) {
			newValue = maskValue;
			if (Rin[r,c] == maskValue) {
				count++;
				# For each cell in the element
				for er = -1 to 1 {
					for ec = -1 to 1 {
						# If the element value isn't null
						if (!IsNull(element[er+2,ec+2])) {

							# If the element is within the bounds of Rin
							if ( ((r + er) >= 1) and ((r + er) <= NumLins(W)) and ((c + ec) >= 1) and ((c + ec) <= NumCols(W)) ) {   

								tempValue = Rin[r + er, c + ec];
								if (tempValue < newValue) {
									newValue = tempValue;
								} # if
							} # if
						} # if
					} # for ec
				} # for er
				if (newValue != maskValue) {
					count--;
					Rtemp[r,c] = newValue;
				} # if
			} # for if
		} # for c
	} # for r
	RasterCopy(Rtemp, Rin);
	print("Erosion Count:", count);
	return count;
} # end MyErosion

######################################################################################################
# This procedure copies the input raster to an output raster and then morphologically
# dilates the masked cells in output raster.  
proc MyDilation(class RASTER Rin, class RASTER Rout, class RASTER Rmask, numeric maskValue, numeric minValue) {

	RasterCopy(Rin, Rout);
	local numeric r, c, er, ec, tempValue, newValue, count=0;
	for r = 1 to NumLins(Rin) {
		for c = 1 to NumCols(Rin) {
			newValue = minValue;
			if (Rmask[r,c] == maskValue) {

				# For each cell in the element
				for er = -1 to 1 {
					for ec = -1 to 1 {
						# If the element value isn't null
						if (!IsNull(element[er+2,ec+2])) {

							# If the element is within the bounds of Rin
							if ( ((r + er) >= 1) and ((r + er) <= NumLins(W)) and ((c + ec) >= 1) and ((c + ec) <= NumCols(W)) ) {   
								tempValue = Rin[r + er, c + ec];
								if (tempValue > newValue) {
									newValue = tempValue;
								} # if
							} # if
						} # if
					} # for ec
				} # for er
				if (newValue != minValue) {
					Rout[r,c] = newValue;
				} # if
			} # if
		} # for c
	} # for r
} # end of MyDilation


######################################################################################################
# This function compares the two input rasters and creates a binary output 
# raster where a value of 1 means that the input raster cell values are different.
proc MyRasterXOR(class RASTER Rin1, class RASTER Rin2, class RASTER Rout) {

	numeric i, j;
	for i = 1 to NumLins(Rin1) {
		for j = 1 to NumCols(Rin1) {

			if (Rin1[i,j] != Rin2[i,j]) {
				Rout[i,j] = 1;
			} else {
				Rout[i,j] = 0;
			}
		}
	}
} # end MyRasterXOR()

######################################################################################################
# This function sets the cells in the raster that symbolizes the new medial lines
# (Rin1) to the average value of the two contour values that created it. 
proc MyMedialLine(class RASTER Rin1, class RASTER Rin2, class RASTER Rmask) {

	numeric i, j;
	for i = 1 to NumLins(Rin1) {
		for j = 1 to NumCols(Rin1) {

			if (Rmask[i,j]) {
				Rin1[i,j] = (Rin2[i,j] + Rin1[i,j]) / 2;
			}
		}
	}
}

######################################################################################################
# This function adds the medial lines that were created to the final output raster
# (also used as the starting raster for the next iteration through the process).  It returns the 
# number of new cell values added to the output raster.
func MyAddToAccumulator(class RASTER Rin, class RASTER Rmask, class RASTER Rout, numeric intContValue) {

	numeric i, j, count = 0;
	for i = 1 to NumLins(Rin) {
		for j = 1 to NumCols(Rin) {

			if ( (Rmask[i,j]) and (Rout[i,j] == intContValue) ) {
				Rout[i,j] = Rin[i,j];
				count++;
			}
		}
	}
	return count;
}

printf ("%s %.2f %.2f %.2f", "Global Mean, Median and Variance values are:", GlobalMean (DEM), GlobalMedian (DEM),GlobalVariance(DEM));