######################################################################################################
#
# 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));
Download