# ------------------------------------------------------------
# OBJECT.sml
# ------------------------------------------------------------
# SET WARNING LEVEL:
$warnings 3
# ------------------------------------------------------------
# COMPUTE-WATERSHED FUNCTION flag STRING ASSIGNMENT:
# See ComputeWatershed SML Function for More Complex Options.
string computeW$ = "Ridge";
# ------------------------------------------------------------
# DEFINE STRING VARIABLES FOR POPUP FUNCTION PROMPTS:
string p$,p1$,p2$,p3$,p4$,p5$,p6$;
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle:
proc writeTitle() begin
printf("DRAFT OBJECT.sml:\n");
printf(" VERSION: May 4, 2006\n");
printf(" PURPOSE: PRODUCE SCENE-OBJECT ");
printf("POLYGONS (SOPs).\n");
printf(" FROM AN INPUT IMAGE RASTER.\n");
printf(" DETAILS: FAQs_by_Jack H");
printf(" AUTHOR: Dr. Jack F. Paris\n");
printf(" CONTACT INFO: [email protected] ");
printf(" 303-775-1195\n");
printf(" ALLOWED USE: ONLY NON-COMMERCIAL\n\n");
end
# ------------------------------------------------------------
# CLEAR CONSOLE WINDOW & ALLOW USER TO REPOSITION IT:
clear();
p1$ = "CONSOLE-WINDOW ADJUSTMENT\n";
p2$ = "* IF NECESSARY, REPOSITION the CONSOLE WINDOW.\n";
p3$ = "* Then, CLICK the OK Button."; p$ = p1$ + p2$ + p3$;
PopupMessage(p$);
# ------------------------------------------------------------
# DECLARE VARIABLES, RASTER OBJECTS, and VECTOR OBJECTS:
class DATABASE db; class DBTABLEINFO table;
class WATERSHED w;
string rObjFileName$,rFilePath$,rFile$,rObj$,rinObj$;
string fNameEP$,fNameW$,fNameSOP$,fNameDum$,desc$;
string fNameEP2$;
string epFile$,epObj$,vNameW$,vOptsW$,wFile$,wObj$;
string rtype$,eptype$,vtopo$,tname$,tdesc$,f1$,f2$;
numeric err,rObj,epObj,linscale,colscale;
numeric optEP;
numeric floorRin,ceilRin,minEP,eLen,cPTLinD,fThin,aIPMin;
numeric lin,col,nlins,ncols,linm1,linp1,colm1,colp1;
numeric linbegin,linend,colbegin,colend,i;
numeric hnRin,nvRin,minRin,maxRin,vR;
numeric r1,r2,r3,r4,r6,r7,r8,r9;
numeric s1,s2,s3,s4,s5,s6,s7,s8,s9,slope;
numeric rflag,rg1,rg2,rg3,rg4,ravg,cIso;
numeric gradx,grady,gradmag;
numeric vEP,vEPold,aminEP,amaxEP,epLin,epFac;
numeric delPar,delParMin,delParMax;
numeric fqsum,cfq,cf10,cf20,cf30,cf40,cf50;
numeric cf60,cf70,cf80,cf90,cfU;
numeric maxEP,eLenp1;
numeric aIPL,aIPR,aIPv,nodeS,nodeE;
numeric numE,numEM,nvpolys,nvlins;
numeric nL,pL,pR,medL,medR,medMin,medMax,medOff,medOff2;
numeric maxmed,ccL,ccR;
numeric ccMin,lenL,lenLMax,wObj;
raster Rin,R,EP,EP2,Dum; vector WATERSHED,SOP;
# ------------------------------------------------------------
# DEFINE proc checkHisto(Band)
# PURPOSE: Ensures that Input Raster Has Unsampled Histogram.
proc checkHisto (raster Band) begin
local numeric sampleInt;
sampleInt = HistogramGetSampleInterval(Band);
if (sampleInt > 1) begin
DeleteHistogram(Band); CreateHistogram(Band,0);
end
else if (sampleInt == -1) then begin
CreateHistogram(Band,0);
end
end
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION:
writeTitle();
# ------------------------------------------------------------
# GET AN INPUT RASTER, Rin:
printf("GOOD CHOICES FOR Rin: Any Tasseled Cap Raster");
printf(" or Any SRFI Raster\n\n");
GetInputRaster(Rin); nlins = NumLins(Rin);
minRin = GlobalMin(Rin); maxRin = GlobalMax(Rin);
rFile$ = GetObjectFileName(Rin);
rObj = GetObjectNumber(Rin);
rinObj$ = GetObjectName(rFile$,rObj);
ncols = NumCols(Rin); rtype$ = RastType(Rin);
checkHisto(Rin); linscale = LinScale(Rin);
colscale = ColScale(Rin); hnRin = HasNull(Rin);
if (hnRin == 1) then nvRin = NullValue(Rin);
# ------------------------------------------------------------
# OPTION: CONTINUE WITH AN EXISTING EP RASTER
# OR: CREATE A NEW EP RASTER:
p1$ = "EP RASTER OPTION:\n";
p2$ = " 1: CREATE a NEW EP Raster\n";
p3$ = " 2: CONTINUE with an EXISTING EP Raster\n";
p4$ = "OPTION ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$;
optEP = PopupNum(p$,1,1,2,0);
# ------------------------------------------------------------
# DEFINE VARIOUS FILE & PATH STRINGS USED IN THIS SCRIPT:
rObjFileName$ = GetObjectFileName(Rin);
rFilePath$ = FileNameGetPath(rObjFileName$);
fNameEP$ = rFilePath$ + "\EP.rvc";
fNameDum$ = rFilePath$ + "\DUMMY.rvc";
fNameW$ = rFilePath$ + "\WATERSHED.rvc";
fNameSOP$ = rFilePath$ + "\SOP.rvc";
# ------------------------------------------------------------
# PRINT INFORMATION ABOUT RASTER Rin TO CONSOLE WINDOW:
printf("INPUT RASTER, Rin, INFORMATION:\n");
printf(" Raster Name: %s\n",rinObj$);
printf(" Data Range: From ");
printf("%6d to %5d\n",minRin,maxRin);
printf(" Line Spacing: %5.2f m ",linscale);
printf("Column Spacing: %5.2f m\n\n",colscale);
# ------------------------------------------------------------
# INITIAL VALUES FOR USER-SPECIFIED CONTROL-PARAMETERS:
floorRin = minRin;
ceilRin = maxRin;
minEP = 1;
eLen = 1;
cPTLinD = 0;
fThin = 0;
aIPMin = 1;
medOff = 0;
# ------------------------------------------------------------
# DEFINITIONS OF CONTROL PARAMETERS:
# floorRin = LOWEST ALLOWED VALUE for RASTER Rin
# ceilRin = HIGHEST ALLOWED VALUE for RASTER Rin
# minEP = THE MINIMUM ALLOWED VALUE FOR RASTER EP
# eLen = THE NUMBER OF CELLS FOR INTERPOLATION
# cPTLinD = CUMULATIVE-PERCENTILE FOR LINE DELETION
# fThin = THE VECTOR-LINE THINNING PARAMETER
# aIPMin = SMALLEST ALLOWED AREA FOR ISLAND POLYGONS
# medOff = ADDED OFFSET FOR SOP MEDIAN VALUES
# ------------------------------------------------------------
# USER INPUTS FOR CONTROL PARAMETERS:
# floorRin:
p1$ = "INPUT-RASTER FLOOR-VALUE:\n";
p2$ = " DEFAULT = MINIMUM VALUE OF RASTER Rin\n";
p3$ = " IF PURPOSE IS TO DEFINE VEG OBJECTS:\n";
p4$ = " THEN USE TC2 AS Rin AND\n";
p5$ = " CHANGE VALUE NEAR ZERO.\n";
p6$ = "VALUE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
floorRin = PopupNum(p$,floorRin,minRin,maxRin,0);
# ceilRin:
p1$ = "INPUT-RASTER CEILING-VALUE:\n";
p2$ = " DEFAULT = MAXIMUM VALUE OF RASTER Rin\n";
p3$ = " IF PURPOSE IS TO DEFINE VEG OBJECTS:\n";
p4$ = " THEN USE TC2 AS Rin AND\n";
p5$ = " CHANGE VALUE NEAR 2000.\n";
p6$ = "VALUE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
ceilRin = PopupNum(p$,ceilRin,minRin,maxRin,0);
# minEP:
p1$ = "MINIMUM ALLOWED EDGE-PROBABILITY VALUE:\n";
p2$ = " RANGE: 1 to 5000\n";
p3$ = "VALUE ENTERED:"; p$ = p1$ + p2$ + p3$;
minEP = PopupNum(p$,minEP,1,5000,0); maxEP = 10000;
# eLen:
p1$ = "EXTRAPOLATION-LINE LENGTH:\n";
p2$ = " RANGE: 1 to 9\n";
p3$ = " LIMIT: MUST BE AN ODD INTEGER\n";
p4$ = " NOTE: ENTER 1 FOR NO EXTRAPOLATION.\n";
p5$ = " PURPOSE: SMOOTHING & INTERPOLATION\n";
p6$ = "VALUE ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ +p6$;
eLen = PopupNum(p$,eLen,1,9,0);
if (eLen > 2) then begin
eLen = 2 * int(eLen / 2) + 1;
end
eLenp1 = eLen + 1;
# cPTLinD:
if (optEP == 1) then begin
p1$ = "CUMULATIVE-PERCENTILE FOR LINE DELETION:\n";
p2$ = " RANGE: 0 to 100 (percent)\n";
p3$ = " NOTE: ENTER 0 FOR NO DELETION.\n";
p3$ = "VALUE ENTERED:"; p$ = p1$ + p2$ + p3$;
cPTLinD = PopupNum(p$,cPTLinD,0,100,0);
end
if (optEP == 2) then cPTLinD = 0;
# fThin:
p1$ = "VECTOR-LINE THINNING-FACTOR:\n";
p2$ = " RANGE: 0.00 to 3.00\n";
p3$ = " NOTE: If < 0.2, NO THINNING.\n";
p4$ = "VALUE ENTERED:"; p$ = p1$ + p2$ + p3$ + p4$;
fThin = PopupNum(p$,fThin,0,3,2);
# aIPMin:
p1$ = "MINIMUM ISLAND-POLYGON AREA:\n";
p2$ = " RANGE: 1 to 1 million sq. m.\n";
p3$ = "VALUE ENTERED:"; p$ = p1$ + p2$ + p3$;
aIPMin = PopupNum(p$,aIPMin,1,1000000,0);
# medOff:
if (cPTLinD > 0) then begin
if (optEP == 1) then begin
p1$ = "OFFSET FOR SOP MEDIAN VALUES:\n";
p2$ = " RANGE: 0 to 10000\n";
p3$ = "VALUE ENTERED:"; p$ = p1$ + p2$ + p3$;
medOff = PopupNum(p$,medOff,0,10000,0);
end
end
# ------------------------------------------------------------
# PRINT USER-SELECTABLE PARAMETERS TO CONSOLE WINDOW:
printf("USER-SELECTABLE CONTROL PARAMETER VALUES:\n");
printf(" floorRin: %7d\n",floorRin);
printf(" ceilRin: %7d\n",ceilRin);
printf(" minEP: %7d\n",minEP);
printf(" eLen: %7d (cells)\n",eLen);
if (optEP == 1) then begin
printf(" cPTLinD: %10.2f (percentile)\n",cPTLinD);
end
printf(" fThin: %10.2f\n",fThin);
printf(" aIPMin: %7d (sq. m.)\n",aIPMin);
if (cPTLinD > 0) then begin
if (optEP == 1) then begin
printf(" medOff: %7d\n",medOff);
end
end
printf("\n");
# ------------------------------------------------------------
# CREATE TEMPORARY RASTER, R, WHICH IS A COPY OF Rin:
CreateTempRaster(R,nlins,ncols,rtype$);
if (hnRin == 1) then IgnoreNull(Rin);
IgnoreNull(R);
for each R begin
vR = Rin; R = vR;
end
CopySubobjects(Rin,R,"GEOREF");
if (hnRin == 1) then begin
SetNull(Rin,nvRin); DeleteHistogram(Rin);
CreateHistogram(Rin,0);
end
rFile$ = GetObjectFileName(R); rObj = GetObjectNumber(R);
rObj$ = GetObjectName(rFile$,rObj);
# ------------------------------------------------------------
# SET RASTER R VALUE TO floorRin OR TO ceilRin:
numE = 0;
for each R begin
vR = R;
if (vR < floorRin) then begin
R = floorRin; numE = numE + 1;
end
if (vR > ceilRin) then begin
R = ceilRin; numE = numE + 1;
end
end
DeleteHistogram(R); CreateHistogram(R,0);
printf("NUMBER OF RASTER R PIXELS CHANGED: %d\n",numE);
# ------------------------------------------------------------
# CREATE A NEW EDGE-PROBABILITIES (EP) RASTER:
desc$ = "Based on " + rinObj$; eptype$ = "16-bit unsigned";
if (optEP == 2) then begin
GetInputRaster(EP,nlins,ncols,eptype$);
printf("SCRIPT CONTINUES WITH AN EXISTING EP RASTER.\n");
end
else begin
CreateRaster(EP,fNameEP$,"EP",desc$,nlins,ncols,eptype$);
printf("SCRIPT STARTS WITH A NEW EP RASTER.\n");
CopySubobjects(R,EP,"GEOREF");
end
epFile$ = GetObjectFileName(EP);
epObj = GetObjectNumber(EP);
epObj$ = GetObjectName(epFile$,epObj);
CreateRaster(Dum,fNameDum$,"Dum","",1,1,"8-bit unsigned");
DeleteObject(fNameDum$,1); PackRVC(fNameDum$);
CreateRaster(EP2,fNameDum$,"EP2","",nlins,ncols,eptype$);
EP2 = minEP;
# ------------------------------------------------------------
# CREATE A NEW SCENE-OBJECT-POLYGONS (SOP) VECTOR:
desc$ = "Scene-Object-Polygons";
CreateVector(SOP,fNameSOP$,"SOP",desc$,"VectorToolkit");
# ------------------------------------------------------------
# SET UP ARRAYS & PARAMETERS FOR EDGE-PROBABILITY PROCESS:
if (eLen == 1) then begin
eLen = 2; eLenp1 = 3;
end
array numeric x[eLenp1]; array numeric y[eLenp1];
linbegin = eLen + 1; colbegin = eLen + 1;
linend = nlins - eLen; colend = ncols - eLen;
for i=1 to eLen begin
x[i] = i - 1;
end
# ------------------------------------------------------------
# SOBEL STAR CONFIGURATION: xx: MARKS RELATIVE POSITIONS
# OF THE EXTENDED "STAR" CELLS USED TO EXTRAPOLATE TO
# THE SOBEL 3 by 3 FILTER ELEMENTS.
#
# THE SOBEL-STAR CONFIGURATION BELOW IS FOR eLen = 3.
#
# xx xx xx
#
# xx xx xx
#
# r1 r2 r3
#
# xx xx r4 cc r6 xx xx
#
# r7 r8 r9
#
# xx xx xx
#
# xx xx xx
#
# ------------------------------------------------------------
# SCAN THE IMAGE TO FIND amaxEP:
printf("SCANNING TO FIND THE MAXIMUM EP VALUE ... ");
amaxEP = 0;
for lin = linbegin to linend begin
for col = colbegin to colend begin
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r1:
for i=1 to eLen y[i] = R[(lin-i),(col-i)];
slope = LinearRegression(x,y,eLen,s1,r1);
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r2:
for i=1 to eLen y[i] = R[(lin-i),col];
slope = LinearRegression(x,y,eLen,s2,r2);
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r3:
for i=1 to eLen y[i] = R[(lin-i),(col+i)];
slope = LinearRegression(x,y,eLen,s3,r3);
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r4:
for i=1 to eLen y[i] = R[lin,(col-i)];
slope = LinearRegression(x,y,eLen,s4,r4);
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r6:
for i=1 to eLen y[i] = R[lin,(col+i)];
slope = LinearRegression(x,y,eLen,s6,r6);
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r7:
for i=1 to eLen y[i] = R[(lin+i),(col-i)];
slope = LinearRegression(x,y,eLen,s7,r7);
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r8:
for i=1 to eLen y[i] = R[(lin+i),col];
slope = LinearRegression(x,y,eLen,s8,r8);
# ESTIMATE VALUE (BY EXTRAPOLATION) AT CELL r9:
for i=1 to eLen y[i] = R[(lin+i),(col+i)];
slope = LinearRegression(x,y,eLen,s9,r9);
# gradx (Spacing-Corrected Sobel 3 by 3 ALGORITHM):
gradx = r3 + 2*r6 + r9 - (r1 + 2*r4 + r7);
gradx = gradx / (8 * colscale);
# grady (Spacing-Corrected Sobel 3 by 3 ALGORITHM):
grady = r7 + 2*r8 + r9 - (r1 + 2*r2 + r3);
grady = grady / (8 * linscale);
# vEP (Absolute Gradiant Vector Magnitude):
vEP = sqrt(gradx*gradx + grady*grady);
if (vEP > amaxEP) then amaxEP = vEP;
end
end
printf("Done\n\n");
epFac = maxEP / amaxEP;
printf("RESULTS OF EP SCAN:\n");
printf(" Before being scaled up by epFac:\n");
printf(" aminEP = %7.2f amaxEP = %10.2f\n",aminEP,amaxEP);
printf(" epFac = %7.2f\n\n",epFac);
# ------------------------------------------------------------
# PROCESS RASTER R TO PRODUCE EDGE-PROBABILITY (EP) VALUES:
printf("PRODUCING FINAL VALUES FOR EP RASTER ... ");
aminEP = 2^16; amaxEP = 0;
for lin = linbegin to linend begin
for col = colbegin to colend begin
for i=1 to eLen y[i] = R[(lin-i),(col-i)];
slope = LinearRegression(x,y,eLen,s1,r1);
for i=1 to eLen y[i] = R[(lin-i),col];
slope = LinearRegression(x,y,eLen,s2,r2);
for i=1 to eLen y[i] = R[(lin-i),(col+i)];
slope = LinearRegression(x,y,eLen,s3,r3);
for i=1 to eLen y[i] = R[lin,(col-i)];
slope = LinearRegression(x,y,eLen,s4,r4);
for i=1 to eLen y[i] = R[lin,(col+i)];
slope = LinearRegression(x,y,eLen,s6,r6);
for i=1 to eLen y[i] = R[(lin+i),(col-i)];
slope = LinearRegression(x,y,eLen,s7,r7);
for i=1 to eLen y[i] = R[(lin+i),col];
slope = LinearRegression(x,y,eLen,s8,r8);
for i=1 to eLen y[i] = R[(lin+i),(col+i)];
slope = LinearRegression(x,y,eLen,s9,r9);
gradx = r3 + 2*r6 + r9 - (r1 + 2*r4 + r7);
gradx = gradx / (8 * colscale);
grady = r7 + 2*r8 + r9 - (r1 + 2*r2 + r3);
grady = grady / (8 * linscale);
gradmag = sqrt(gradx*gradx + grady*grady);
vEP = ceil(epFac * gradmag);
if (vEP < aminEP) then aminEP = vEP;
if (vEP > amaxEP) then amaxEP = vEP;
vEP = Bound(vEP,minEP,maxEP);
if (optEP == 2) then begin
vEPold = EP[lin,col];
if (vEP > vEPold) then begin
EP[lin,col] = vEP;
end
end
else EP[lin,col] = vEP;
end
end
# ------------------------------------------------------------
# FIND AND REPLACE ISOLATED EP CELLS THAT HAVE VALUES = minEP.
# REPLACEMENT VALUE IS THE MAXIMUM OF THE FOUR AVERAGES
# OF OPPOSING CELLS (HORIZONTAL, VERTICAL & THE TWO DIAGONALS):
cIso = 0;
for lin = linbegin to linend begin
for col = colbegin to colend begin
vEP = EP[lin,col];
if (vEP == minEP) then begin
rflag = 0; rg1 = 0; rg2 = 0; rg3 = 0; rg4 = 0;
linm1 = lin - 1; linp1 = lin + 1;
colm1 = col - 1; colp1 = colp1;
r1 = EP[linm1,colm1];
r9 = EP[linp1,colp1];
if (r1 > minEP and r9 > minEP) then begin
rflag = 1; rg1 = (r1 + r9) / 2;
end
r2 = EP[linm1,col];
r8 = EP[linp1,col];
if (r2 > minEP and r8 > minEP) then begin
rflag = 1; rg2 = (r2 + r8) / 2;
end
r3 = EP[linm1,colp1];
r7 = EP[linp1,colm1];
if (r3 > minEP and r7 > minEP) then begin
rflag = 1; rg3 = (r3 + r7) / 2;
end
r4 = EP[lin,colm1];
r6 = EP[lin,colp1];
if (r4 > minEP and r6 > minEP) then begin
rflag = 1; rg4 = (r4 + r6) / 2;
end
if (rflag == 1) then begin
ravg = SetMax(rg1,rg2,rg3,rg4);
EP[lin,col] = ceil(ravg);
cIso = cIso + 1;
end
end
end
end
if (eLen == 2) then eLen = 1;
for lin = linbegin to linend begin
for col = colbegin to colend begin
EP2[lin,col] = FocalMax(EP,eLen,eLen,lin,col,1);
end
end
# TRANSFER FILTERED EP2 VALUES TO EP RASTER:
EP = EP2;
printf("Done\n");
printf("RASTER EP INFORMATION:\n");
printf(" EP RANGE: From %4d to %7d\n",minEP,maxEP);
printf(" Isolated EP Cells (Changed): %d\n\n",cIso);
IgnoreNull(EP); DeleteHistogram(EP); CreateHistogram(EP,0);
# ------------------------------------------------------------
# ASSIGN HANDLE (w) to EP RASTER:
# NOTE: EP is treated as a "DEM" so that the TNTmips
# WATERSHED processes can be used to produce
# Scene-Object Polygons (i.e., "WATERSHEDS").
w = WatershedInit(epFile$,epObj$);
# ------------------------------------------------------------
# PROCESS RASTER EP TO MAKE SCENE-OBJECT POLYGONS (SOPs):
printf("PROCESSING RASTER EP TO PRODUCE POLYGONS:\n");
printf(" Making Scene-Object Polygons (SOPs) .... ");
WatershedCompute(w,computeW$); printf("Done\n");
printf(" Saving SOPs to WATERSHED Vector ........ ");
WatershedGetObject(w,"VectorWatershed",wFile$,wObj$);
wObj = ObjectNumber(wFile$,wObj$,"Vector");
CopyObject(wFile$,wObj,fNameDum$);
WatershedClose(w); CreateHistogram(EP,0);
CreatePyramid(EP,0); CloseRaster(EP);
CloseRaster(EP2); printf("Done\n\n");
# ------------------------------------------------------------
# RE-OPEN VECTOR OBJECT, WATERSHED, FOR FURTHER PROCESSING:
vNameW$ = "WATERSHED"; vOptsW$ = "VectorToolkit";
OpenVector(WATERSHED,fNameDum$,vNameW$,vOptsW$);
numEM = WATERSHED.$Info.NumPoints;
# ------------------------------------------------------------
# DELETE POUR POINTS FROM WATERSHED VECTOR OBJECT:
array numeric ptlist[numEM + 1];
for numE=1 to numEM begin
ptlist[numE] = numE;
end
if (numEM > 0) then begin
VectorDeletePoints(WATERSHED,ptlist,numEM);
end
# ------------------------------------------------------------
# COPY WATERSHED VECTOR ELEMENTS TO VECTOR OBJECT SOP:
nvlins = WATERSHED.$Info.NumLines;
nvpolys = WATERSHED.$Info.NumPolys;
VectorCopyElements(WATERSHED,SOP);
CloseVector(WATERSHED);
# ------------------------------------------------------------
# DELETE LARGE POLYGON AROUND THE OUTSIDE EDGE OF RASTER EP:
numeric polyNum, numPolyMaxArea, area, maxArea;
for polyNum=1 to nvpolys begin
area = SOP.poly[polyNum].POLYSTATS.Area;
if (area > maxArea) begin
maxArea = area;
numPolyMaxArea = polyNum;
end
end
if (SOP.Poly[numPolyMaxArea].Internal.MinX == 0 and SOP.Poly[numPolyMaxArea].Internal.MinY == 0) then
VectorDeletePoly(SOP, numPolyMaxArea);
# ------------------------------------------------------------
# BUILD POLYGONS FROM MODIFIED LINES:
printf("BUILDING INITIAL SOPs IN VECTOR SOP ....... ");
VectorValidate(SOP);
printf("Done\n");
# ------------------------------------------------------------
# MAKE A REPORT CONCERNING VECTOR SOP's ELEMENTS:
nvpolys = SOP.$Info.NumPolys;
nvlins = SOP.$Info.NumLines;
printf(" Polygons: %d",nvpolys);
printf(" Lines: %d\n\n",nvlins);
printf("REMAINING PROCESSING STEPS");
printf(" STATUS\n");
# ------------------------------------------------------------
# COMPUTE RASTER PROPERTIES (FROM R) FOR R_Props TABLE:
tname$ = "R_Props"; tdesc$ = "R Attributes";
f1$ = "Proportionally"; f2$ = "NULL";
array numeric delLlist[nvlins];
# ------------------------------------------------------------
if (optEP == 1) then begin
# START PROCESS TO COMBINE & CLEAN UP SIMILAR SOPs:
printf(" CREATING R_Props DATABASE TABLE ........ ");
ComputeRasterProperties(Rin,SOP,tname$,tdesc$,f1$,f2$);
printf("Done\n ");
if (cPTLinD > 0) then begin
# PROCESS TO DELETE LINES BETWEEN SIMILAR SOPs:
# DETERMINE RANGE OF DELETION PARAMETER (delPar):
printf("DETERMINING RANGES ..................... ");
medMin = 100000; medMax = -100000;
for nL=1 to nvlins begin
pL = SOP.line[nL].Internal.LeftPoly;
pR = SOP.line[nL].Internal.RightPoly;
if (pL > 0 and pR > 0) then begin
ccL = SOP.poly[pL].R_Props.CellCount;
ccR = SOP.poly[pR].R_Props.CellCount;
ccMin = SetMin(ccL,ccR);
if (ccMin > 5) then begin
medL = SOP.poly[pL].R_Props.Median;
medR = SOP.poly[pR].R_Props.Median;
if (medL < medMin) then medMin = medL;
if (medL > medMax) then medMax = medL;
if (medR < medMin) then medMin = medR;
if (medR > medMax) then medMax = medR;
end
end
end
medOff2 = medOff - medMin + 10;
delParMin = 100000; delParMax = -100000;
for nL = 1 to nvlins begin
pL = SOP.line[nL].Internal.LeftPoly;
pR = SOP.line[nL].Internal.RightPoly;
if (pL > 0 and pR > 0) then begin
ccL = SOP.poly[pL].R_Props.CellCount;
ccR = SOP.poly[pR].R_Props.CellCount;
ccMin = SetMin(ccL,ccR);
if (ccMin > 5) then begin
medL = SOP.poly[pL].R_Props.Median;
medR = SOP.poly[pR].R_Props.Median;
medL = medL + medOff2; medR = medR + medOff2;
delPar = abs((medL - medR) / (medL + medR));
delPar = round(10000 * delPar);
if (delPar < delParMin) then delParMin = delPar;
if (delPar > delParMax) then delParMax = delPar;
end
end
end
printf("Done\n");
printf(" Polygon Rin Median Range: %6d",medMin);
printf(" to %6d\n",medMax);
printf(" delPar Range: %7d",delParMin);
printf(" to %7d\n",delParMax);
# CALCULATE CUMULATIVE DISTRIBUTION OF delPar VALUES:
array numeric fq[delParMax + 1];
# CALCULATE RAW FREQUENCY DISTRIBUTION (fq[i]):
for nL = 1 to nvlins begin
pL = SOP.line[nL].Internal.LeftPoly;
pR = SOP.line[nL].Internal.RightPoly;
if (pL > 0 and pR > 0) then begin
ccL = SOP.poly[pL].R_Props.CellCount;
ccR = SOP.poly[pR].R_Props.CellCount;
ccMin = SetMin(ccL,ccR);
if (ccMin > 5) then begin
medL = SOP.poly[pL].R_Props.Median;
medR = SOP.poly[pR].R_Props.Median;
medL = medL + medOff2; medR = medR + medOff2;
delPar = abs((medL - medR) / (medL + medR));
delPar = round(10000 * delPar);
fq[delPar] = fq[delPar] + 1;
end
end
end
fqsum = 0;
# CALCULATE SUM OF RAW FREQUENCY VALUES:
for i=0 to delParMax begin
fqsum = fqsum + fq[i];
end
# CALCULATE CUMULATIVE FREQUENCY DISTRIBUTION:
cfq = 0; cf10 = 0; cf20 = 0; cf30 = 0; cf40 = 0;
cf50 = 0; cf60 = 0; cf70 = 0; cf80 = 0; cfU = 0;
for i=0 to delParMax begin
cfq = cfq + fq[i] / fqsum;
fq[i] = 100 * cfq;
if (cfU == 0) then begin
if (fq[i] >= cPTLinD) then cfU = i;
end
if (cf10 == 0) then begin
if (fq[i] >= 10) then cf10 = i;
end
if (cf20 == 0) then begin
if (fq[i] >= 20) then cf20 = i;
end
if (cf30 == 0) then begin
if (fq[i] >= 30) then cf30 = i;
end
if (cf40 == 0) then begin
if (fq[i] >= 40) then cf40 = i;
end
if (cf50 == 0) then begin
if (fq[i] >= 50) then cf50 = i;
end
if (cf60 == 0) then begin
if (fq[i] >= 60) then cf60 = i;
end
if (cf70 == 0) then begin
if (fq[i] >= 70) then cf70 = i;
end
if (cf80 == 0) then begin
if (fq[i] >= 80) then cf80 = i;
end
if (cf90 == 0) then begin
if (fq[i] >= 90) then cf90 = i;
end
end
# PRINTOUT CUMULATIVE FREQUENCY DISTRIBUTION DATA:
printf(" delPar at %4.1f Percentile:",cPTLinD);
printf(" %d \n",cfU);
printf(" CUMULATIVE DISTRIBUTION OF delPar:\n");
printf(" cfq delPar\n");
printf(" 10%9d\n",cf10);
printf(" 20%9d\n",cf20);
printf(" 30%9d\n",cf30);
printf(" 40%9d\n",cf40);
printf(" 50%9d\n",cf50);
printf(" 60%9d\n",cf60);
printf(" 70%9d\n",cf70);
printf(" 80%9d\n",cf80);
printf(" 90%9d\n",cf90);
# BEGIN PROCESS TO DELETE WEAK POLYLINES:
printf(" FINDING WEAK EDGES ..................... ");
for nL=1 to nvlins begin
delLlist[nL] = 0;
end
numEM = 0;
for nL = 1 to nvlins begin
pL = SOP.line[nL].Internal.LeftPoly;
pR = SOP.line[nL].Internal.RightPoly;
if (pL > 0 and pR > 0) then begin
ccL = SOP.poly[pL].R_Props.CellCount;
ccR = SOP.poly[pR].R_Props.CellCount;
ccMin = SetMin(ccL,ccR);
if (ccMin > 5) then begin
medL = SOP.poly[pL].R_Props.Median;
medR = SOP.poly[pR].R_Props.Median;
medL = medL + medOff2; medR = medR + medOff2;
delPar = abs((medL - medR) / (medL + medR));
delPar = round(10000 * delPar);
if (delPar < cfU) then begin
numEM = numEM + 1;
delLlist[numEM] = nL;
end
end
end
end
printf("Done\n ");
printf("COMBINING SOPs THAT ARE ALIKE .......... ");
if (numEM > 0) then begin
VectorDeleteLines(SOP,delLlist,numEM);
end
printf("Done\n ");
# (RE)BUILD SOPs:
printf("(RE)BUILDING SOPs ...................... ");
VectorValidate(SOP);
printf("Done\n ");
end
# BEGIN PROCESS TO DELETE DANGLING LINES:
lenLMax = 0; nvlins = SOP.$Info.NumLines;
for nL=1 to nvlins begin
lenL = SOP.line[nL].LINESTATS.Length;
if (lenL > lenLMax) then lenLMax = lenL;
end
lenLMax = lenLMax + 1;
printf("DELETING DANGLING LINES IN SOP ......... ");
VectorDeleteDangleLines(SOP,lenLMax);
printf("Done\n ");
# RE-BUILD SOPs:
printf("(RE)BUILDING SOPs ...................... ");
VectorValidate(SOP);
printf("Done\n");
end
# ------------------------------------------------------------
# DELETE LINES ASSOCIATED WITH SMALL ISLAND POLYGONS in SOP:
printf(" FINDING ISLAND-POLYGONS TO DELETE ...... ");
nvlins = SOP.$Info.NumLines;
# RE-SET delLlist ARRAY VALUES TO ZERO:
for numE=1 to nvlins begin
delLlist[numE] = 0;
end
numEM = 0;
# LINES ASSOCIATED WITH ISLAND POLYGONS USE SAME NODE NUMBER
# FOR START NODE AND ENDING NODE:
for nL=1 to nvlins begin
nodeS = SOP.line[nL].Internal.StartNode;
nodeE = SOP.line[nL].Internal.EndNode;
# FIND AREA OF SMALLER RELATED POLYGON:
if (nodeS == nodeE) then begin
pL = SOP.line[nL].Internal.LeftPoly;
pR = SOP.line[nL].Internal.RightPoly;
aIPL = SOP.poly[pL].POLYSTATS.Area;
aIPR = SOP.poly[pR].POLYSTATS.Area;
aIPv = SetMin(aIPL,aIPR);
# TEST aIPv (AREA) AGAINST aIPMin:
if (aIPv < aIPMin) then begin
numEM = numEM + 1;
delLlist[numEM] = nL;
end
end
end
# END OF PROCESS TO DEFINE & DELETE WEAK EDGES
# ------------------------------------------------------------
printf("Done\n ");
# PERFORM THE LINE-DELETION PROCESS (FOR SMALL ISLANDS):
if (numEM > 0) then begin
printf("DELETING %5d ISLAND-POLYGONS ......... ",numEM);
if (numEM > 0) then begin
VectorDeleteLines(SOP,delLlist,numEM);
end
printf("Done\n ");
printf("(RE)BUILDING SOPs ...................... ");
VectorValidate(SOP);
printf("Done\n ");
end
if (numEM == 0) then begin
printf(" No Island Polygons Were Deleted.\n ");
end
# ------------------------------------------------------------
if (optEP == 1) then begin
# DELETE TABLE NAMED R_Props:
db = OpenVectorPolyDatabase(SOP);
table = DatabaseGetTableInfo(db,tname$);
printf("DELETING R_Props DATABASE TABLE ........ ");
err = TableDelete(table);
printf("Done\n ");
end
# IF optEP == 2, THEN SUBJECT TABLE DOES NOT EXIST.
# ------------------------------------------------------------
if (fThin >= 0.2) then begin
# RUN LINE THINNING PROCESS:
printf("THINNING VECTOR SOP LINES .............. ");
VectorThinLines(SOP,"Douglas-Peucker",fThin);
printf("Done\n ");
VectorRemoveExcessNodes(SOP);
printf("(RE)BUILDING SOPs ...................... ");
VectorValidate(SOP);
printf("Done\n ");
end
# ------------------------------------------------------------
# FINAL DELETION OF DANGLING LINES IN VECTOR SOP:
lenLMax = 0; nvlins = SOP.$Info.NumLines;
for nL=1 to nvlins begin
lenL = SOP.line[nL].LINESTATS.Length;
if (lenL > lenLMax) then lenLMax = lenL;
end
lenLMax = lenLMax + 1;
printf("DELETING DANGLING LINES IN SOP ......... ");
VectorDeleteDangleLines(SOP,lenLMax);
printf("Done\n ");
# ------------------------------------------------------------
# RE-SET Rin TO INPUT CONDITIONS & CLOSE IT:
if (hnRin == 1) then begin
SetNull(Rin,nvRin);
DeleteHistogram(Rin); CreateHistogram(Rin,0);
end
CloseRaster(Rin);
# ------------------------------------------------------------
# BUILDING FINAL SOPs:
printf("(RE)BUILDING FINAL SOPs ................ ");
VectorValidate(SOP);
printf("Done\n\n");
# ------------------------------------------------------------
# MAKE REPORT OF VECTOR ELEMENTS IN FINAL SOP VECTOR OBJECT:
nvlins = SOP.$Info.NumLines;
nvpolys = SOP.$Info.NumPolys;
printf("FINAL VECTOR SOPs & EDGE LINES:\n");
printf(" Polygons: %d",nvpolys);
printf(" Lines: %d\n\n",nvlins);
# ------------------------------------------------------------
# CLOSE VECTOR OBJECT SOP:
CloseVector(SOP);
# ------------------------------------------------------------
# DELETE TEMPORARY RASTER R:
DeleteTempRaster(R);
# ------------------------------------------------------------
printf("TO SAVE THE CONSOLE WINDOW TEXT AS A REPORT:\n");
printf(" 1. RIGHT-CLICK IN THE CONSOLE WINDOW.\n");
printf(" 2. SELECT THE Save As... OPTION.\n");
printf(" 3. NAVIGATE TO THE DESIRED LOCATION.\n");
printf(" 4. PROVIDE A REPORT NAME (or OVERWRITE).\n");
printf(" 5. CLICK OK.\n\n");
printf("*** TO CONTINUE, YOU MUST EXIT SML EDITOR ***");
# ------------------------------------------------------------
# DELETE DUMMY.rvc FILE:
DeleteFile(fNameDum$);
Download