# ------------------------------------------------------------
# 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: jparis37@msn.com ");
     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$    = "Proportionality"; 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$);