# ------------------------------------------------------------
# GRUVI.sml
# ------------------------------------------------------------
# SET WARNING LEVEL:
  $warnings 3;
# ------------------------------------------------------------
# DECLARE SEVERAL VARIABLES RELATED TO GRUVI.sml:
  numeric srfiXorg,srfiYorg,slope,bnp,bnpxy,tsrfi,cvi;
  numeric bnpop,bnpi,bnpf,dbnpf;
  numeric srfiXb,srfiYb,srfiXv,srfiYv;
  numeric srfiX,srfiY;
  numeric tsrfiX1,tsrfiY1,tsrfiX2,tsrfiY2,tsrfiX3,tsrfiY3;
  numeric tsrfiX4,tsrfiY4,tsrfiX5,tsrfiY5,tsrfiX6,tsrfiY6;
  numeric nlins,ncols,isnull,i;
  numeric rang,rangdeg,sinrang,cosrang;
  numeric grubif,grubifmin,grubi,grumax,grunull;
  numeric gmean,gstd,gcv1,gcv2,gflag;
  numeric gruvif,gruvi;
  numeric gruvirange,gruvidveg,gfac;
  numeric hnmGB,hnmGF,hnGB,hnGF,nvGB,nvGF,pGB,pGF;
  numeric c1GB,c2GB,c3GB,c1GF,c2GF,c3GF;
  numeric pVI;
  raster SRFIX,SRFIY,GRUVI,GRUBI,GRUFI;
  raster SRFIXT,SRFIYT,GRUVIT;
  string gtype$,rtype$,fore$,back$;
# ------------------------------------------------------------
# DECLARE VARIABLES RELATED TO CONTRAST LOOKUP PROCESS: 
  class CONTRAST smlContrast;
  string clut1$,clut2$,clut3$;
  array numeric v[131072];
  clut1$  = "exp"; clut2$  = "Exponential";
  clut3$  = "Exponential contrast table";
  smlContrast.OutputLowerLimit =   1;
  smlContrast.OutputUpperLimit = 255;
# ------------------------------------------------------------
# DECLARE SEVERAL VARIABLES FOR GENERAL CHARACTER STRINGS:
  string t$,p$,p1$,p2$,p3$,p4$,p5$,p6$,p7$,p8$,p9$;
  string p10$,p11$,p12$,p13$,p14$,p15$,p16$,p17$,p18$,p19$;
  string p20$,p21$;
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle:
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
  proc writeTitle(pVIloc) begin
     numeric pVIloc;
     printf("GRUVI.sml:\n");
     printf("         VERSION: November 15, 2005\n");
     if (pVIloc == 1) then begin
        printf("         PURPOSE: PRODUCE GRUVI & GRUBI ");
        printf("RASTERS FROM A PAIR OF\n");
        printf("                  SRFIX & SRFIY RASTERS.\n");
     end
     else begin
        printf("         PURPOSE: PRODUCE GRUFI & GRUBI ");
        printf("RASTERS FROM A PAIR OF\n");
        printf("                  SRFIX & SRFIY RASTERS,\n");
        printf("                  OR FROM A PAIR OF TCX");
        printf(" AND TCY RASTERS (FROM TASCAP.SML).\n");
     end
     printf("     NEW FEATURE: MAKES AUTOMATIC CONTRAST ");
     printf("TABLES.\n");
     printf("         DETAILS: FAQs_by_Jack A & E\n");
     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 & REQUEST REPOSITIONING:
  clear();
  p1$ = "CONSOLE-WINDOW ADJUSTMENT\n";
  p2$ = "* REPOSITION the CONSOLE WINDOW.\n";
  p3$ = "* Then, CLICK OK.";
  p$  = p1$ + p2$ + p3$;
  PopupMessage(p$);
# ------------------------------------------------------------
# SELECT USE TYPE:
  p1$ = "GRUVI.sml can be Used in TWO Ways:\n";
  p2$ = "1: Produce GRUVI & GRUBI (Veg. & Bright.)\n";
  p3$ = "2: Produce GRUFI & GRUBI (Foreground & Bright.)\n";
  p4$ = "   DEFAULT: 1 (Veg. & Bright.)\n";
  p5$ = "   WARNING: If Using TC Rasters, SELECT WAY 2!\n";
  p6$ = "WAY SELECTED (1 or 2):";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
  pVI = PopupNum(p$,1,1,2,0);
# ------------------------------------------------------------
# func GlobalPtile(Rl,ptilel):
  func GlobalPtile(raster Rl, numeric ptilel,hnl,nvl) begin
     local numeric funcval,rl,roffl;
     local numeric nlinsl,ncolsl,linl,coll,il;
     local numeric countl,suml,flagl,cfql;
     roffl = -65535; countl = 0;
     if (hnl == 1) then begin 
        for each Rl begin
           rl = Rl[linl,coll];
           if (rl <> nvl) then countl = countl + 1;
        end
     end
     if (hnl == 0) then begin
        nlinsl = NumLins(Rl); ncolsl = NumCols(Rl);
        countl = nlinsl * ncolsl;
     end
     for il = 0 to 131071 v[il] = 0;
     for each Rl[linl,coll] begin
        rl = Rl[linl,coll];
        if (hnl == 1) then begin
           if (rl <> nvl) then begin
              il = rl - roffl; v[il] = v[il] + 1;
           end
        end
        if (hnl == 0) then begin
           il = rl - roffl; v[il] = v[il] + 1;
        end
     end
     suml = v[0]; flagl = 0;
     for il=1 to 131071 begin
        if (flagl == 0) then begin
           suml = suml + v[il]; cfql = suml * 100 / countl;
           if (cfql >= ptilel) then begin
              flagl = 1; funcval = il + roffl;
           end
        end
     end
     return (funcval);
  end
# ------------------------------------------------------------
# DEFINE func expPower(dnLl,dnMl,dnHl):
  func expPower(dnLl,dnMl,dnHl) begin
     local numeric delLHl,delLMl,sr1l,sr2l,exppl;
     sr1l = 0.3; sr2l = 2 * (1 - sr1l);
     delLHl = dnHl - dnLl; delLMl = dnMl - dnLl;
     exppl = sr1l + sr2l * delLMl / delLHl;
     return (exppl);
  end
# ------------------------------------------------------------
# DEFINE proc checkHisto(Band)
# PURPOSE: Check that Input Raster has a Full (Unsampled)
#          Histogram.  If Not, Then Recompute 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
# ------------------------------------------------------------
# ASSIGN VALUES TO GRUVI DATATYPE, SCALE & RANGE PARAMETERS:
  gtype$ = "16-bit signed";
  grumax = 3000; gruvirange = 1000; grunull = -grumax - 1;
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION:
  writeTitle(pVI);
# ------------------------------------------------------------
# ASK USER TO NAME FOREGROUND & BACKGROUND MATERIALS:
  if (pVI == 1) then begin
     printf("GRUVI RELATES TO FOREGROUND VEGETATION");
     printf(" CONTRASTED AGAINST\nA BACKGROUND OF ");
     printf("BARE SOILS, RESIDUE, OR WATER.\n\n");
  end
  else begin
     t$  = "IDENTIFY FOREGROUND MATERIALS";
     p$ = "FOREGROUND-MATERIALS NAME?";
     fore$ = PopupString(p$,"ENTER NAME",t$);
     printf("  GRUFI RELATES TO: ");
     printf("%s\n",fore$);
     t$  = "IDENTIFY BACKGROUND MATERIALS";
     p$ = "BACKGROUND-MATERIALS NAME?";
     back$ = PopupString(p$,"ENTER NAME",t$);
     printf("CONTRASTED AGAINST: ");
     printf("%s\n\n",back$);
  end
# ------------------------------------------------------------
# ASK USER TO SELECT TRANSFORMATION OPTION:
  if (pVI == 2) then begin
     tsrfi = 1;
  end
  else begin
     p1$ = "SRFI TRANSFORMATION OPTION:\n";
     p2$ = "   For Classic NDVI, SAVI, or PVI: CLICK  NO\n";
     p3$ = "   For Transformed NDVI or  TSAVI: CLICK YES\n";
     p4$ = "             For CUSTOMIZED GRUVI: CLICK YES\n";
     p5$ = "CLICK YES or NO:";
     p$ = p1$ + p2$ + p3$ + p4$ + p5$;
     tsrfi = PopupYesNo(p$,1);
  end
# ------------------------------------------------------------
# ASK USER TO DEFINE FOREGROUND-MATERIALS POINT IN 2-SPACE:
# srfiXv = SRFI VALUE for FOREGROUND MATERIALS.
  if (pVI == 1) then begin
     p1$ = "X-AXIS DENSE-VEGETATION POINT:\n";
     p2$ = "   DEFAULT: 300\n";
     p3$ = "VALUE ENTERED:";
     p$ = p1$ + p2$ + p3$;
     srfiXv = PopupNum(p$,300,-32000,32000,0);
#    srfiYv = SRFIY VALUE for FOREGROUND MATERIALS.
     p1$ = "Y-AXIS DENSE-VEGETATION POINT:\n";
     p2$ = "   DEFAULT: 6000\n";
     p3$ = "VALUE ENTERED:";
     p$ = p1$ + p2$ + p3$;
     srfiYv = PopupNum(p$,6000,-32000,32000,0);
  end
  else begin
     p1$ = "X-AXIS FOREGROUND-MATERIALS POINT:\n";
     p2$ = "   DEFAULT: 300\n";
     p3$ = "VALUE ENTERED:";
     p$ = p1$ + p2$ + p3$;
     srfiXv = PopupNum(p$,300,-32000,32000,0);
#    srfiYv = SRFIY VALUE for FOREGROUND MATERIALS.
     p1$ = "Y-AXIS FOREGROUND-MATERIALS POINT:\n";
     p2$ = "   DEFAULT: 6000\n";
     p3$ = "VALUE ENTERED:";
     p$ = p1$ + p2$ + p3$;
     srfiYv = PopupNum(p$,6000,-32000,32000,0);
  end 
  if (tsrfi == 1) then begin
#    ASK USER TO DEFINE SLOPE OF LBM:
#    LBM is the LINE OF BACKGROUND MATERIALS.
     p1$ = "SLOPE OF LINE-OF-BACKGROUND-MATERIALS:\n";
     p2$ = "VALUE ENTERED:";
     p$ = p1$ + p2$;
     slope = PopupNum(p$,1.086,-100,100,4);
#    ASK USER TO DEFINE BACKGROUND-MATERIALS POINT IN 2-SPACE:
#    srfiXb = SRFIX VALUE for BACKGROUND MATERIALS.
     p1$ = "X-AXIS BACKGROUND-MATERIALS POINT:\n";
     p2$ = "VALUE ENTERED:\n";
     p$ = p1$ + p2$;
     srfiXb = PopupNum(p$,1711,-32000,32000,0);
#    srfiYb = SRFI VALUE for BACKGROUND MATERIALS.
     p1$ = "Y-AXIS BACKGROUND-MATERIALS POINT:\n";
     p2$ = "VALUE ENTERED:";
     p$ = p1$ + p2$;
     srfiYb = PopupNum(p$,2132,-32000,32000,0);
     srfiXorg = 0;
#    IF pVI == 2, ASK USER TO DEFINE NEW X ORIGIN IN 2-SPACE:
#    Xorg = SRFIX VALUE for NEW ORIGIN.
#    Corresponding Yorg is CALCULATED -- not Specified.
     if (pVI == 2) then begin
        p1$ = "LINE-OF-BACKGROUND-MATERIALS X ORIGIN:\n";
        p2$ = "VALUE ENTERED:";
        p$ = p1$ + p2$;
        srfiXorg = PopupNum(p$,0,-32000,32000,0);
     end
#    CALCULATED srfiYorg VALUE:  
     srfiYorg = round(srfiYb - slope * (srfiXb - srfiXorg));
  end
  else begin
     srfiXorg = 0; srfiYorg = 0; slope = 1;
  end
# ------------------------------------------------------------
# ASK USER TO OPT FOR OPTIMIZING THE bnp PARAMETER:
# IF NOT OPTIMIZED, THEN USER MUST SPECIFY THE VALUE OF bnp.
  p1$ = "BACKGROUND NOISE PARAMETER (bnp) OPTIMIZATION:\n";
  p2$ = "Are You USING a TEST AREA for bnp OPTIMIZATION?\n";
  p$ = p1$ + p2$;
  bnpop = PopupYesNo(p$,0);
  if (bnpop == 1) then begin
#    START OPTIMIZATION PROCESS LEADING TO bnp:
     gcv1 = 1000; gflag = 0;
     GetInputRaster(SRFIXT);
     nlins = NumLins(SRFIXT); ncols = NumCols(SRFIXT);
     rtype$ = RastType(SRFIXT);
     checkHisto(SRFIXT);
     GetInputRaster(SRFIYT,nlins,ncols,rtype$);
     checkHisto(SRFIYT);
     rang = -atan(slope);
     rangdeg = rang * deg; rangdeg = rangdeg + 45;
     printf("   CCW ROTATION ANGLE: %6.2f deg\n\n",rangdeg);    
     sinrang = sind(rangdeg); cosrang = cosd(rangdeg);
#    TRANSLATION TO A NEW ORIGIN:
     tsrfiX1 = srfiXv - srfiXorg;
     tsrfiY1 = srfiYv - srfiYorg;
#    SCALE-PRESERVING ROTATION:
     tsrfiX2 = tsrfiX1 * cosrang - tsrfiY1 * sinrang;
     tsrfiY2 = tsrfiX1 * sinrang + tsrfiY1 * cosrang;
     dbnpf = 0.02; bnpf = 0.00;
     CreateTempRaster(GRUVIT,nlins,ncols,gtype$,0);
     IgnoreNull(GRUVIT);
     for i = 0 to 1000 step 1 begin
        if (gflag == 0) then begin
           bnpf = bnpf + dbnpf;
           bnpi = bnpf * 5000;
           tsrfiX3 = tsrfiX2 + bnpi;
           tsrfiY3 = tsrfiY2 + bnpi;
           gruvidveg = (tsrfiY3 - tsrfiX3);
           gruvidveg = gruvidveg / (tsrfiY3 + tsrfiX3); 
           gfac = gruvirange / gruvidveg;
           for each GRUVIT begin
              GRUVIT = grunull;
              isnull = IsNull(SRFIXT);
              if (isnull == 0) then begin
                 srfiX = SRFIXT;
                 srfiY = SRFIYT;
                 tsrfiX4 = srfiX - srfiXorg;
                 tsrfiY4 = srfiY - srfiYorg;
                 tsrfiX5 = tsrfiX4 * cosrang - tsrfiY4 * sinrang;
                 tsrfiY5 = tsrfiX4 * sinrang + tsrfiY4 * cosrang;
                 tsrfiX6 = tsrfiX5 + bnpi;
                 tsrfiY6 = tsrfiY5 + bnpi;
                 gruvif = tsrfiY6 - tsrfiX6;
                 gruvif = gfac * gruvif / (tsrfiY6 + tsrfiX6);
                 gruvi = round(gruvif);
                 if (gruvi < -grumax) then gruvi = -grumax;
                 if (gruvi > grumax) then gruvi = grumax;
                 GRUVIT = gruvi;
              end
           end
           SetNull(GRUVIT,grunull);
           gmean = GlobalMean(GRUVIT);
           gmean = gmean - 1000;
           gstd = GlobalSD(GRUVIT);
           gcv2 = gstd * 100 / gmean;
           if (gcv2 < gcv1) then begin
              gcv1 = gcv2;
           end
           else begin
              bnp = bnpf - dbnpf;
              gflag = 1;
           end
           IgnoreNull(GRUVIT);
        end
     end
     CloseRaster(GRUVIT);
     CloseRaster(SRFIXT); CloseRaster(SRFIYT);
     printf("   OPTIMIZED bnp VALUE: %5.2f\n\n",bnp);
  end
  else bnp = 0.5;
# ASK USER TO SPECIFY THE VALUE OF bnp:
  if (tsrfi == 1) then begin
     p1$ = "BACKGROUND-NOISE PARAMETER (bnp) FOR\n";
     p2$ = "TRANSFORMED INDICES ONLY:\n";
     p3$ = "   RANGE: 0.02 to 20.00\n";
     p4$ = "   Classic    TSAVI: bnp  =  0.15\n";
     p5$ = "   Classic     WDVI: bnp  = 20.00\n";
     p6$ = "   Customized INDEX: 0.00 < bnp < 20.00\n";
     p7$ = "ENTER bnp VALUE:\n";
     p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$;
     bnp = PopupNum(p$,bnp,0.02,20.00,2);
  end
  else begin
     p1$ = "BACKGROUND-NOISE PARAMETER (BNP) FOR\n";
     p2$ = "NON-TRANSFORMED INDICES ONLY:\n";
     p3$ = "   RANGE: 0.0 to 20.00\n";
     p4$ = "   Classic     NDVI: bnp  =  0.00\n";
     p5$ = "   Classic     SAVI: bnp  =  0.50\n";
     p6$ = "   Classic      PVI: bnp  = 20.00\n";
     p7$ = "ENTER bnp VALUE:\n";
     p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$;
     bnp = PopupNum(p$,bnp,0.0,20.00,2);
  end
  printf("   LINE OF BACKGROUND MATERIALS:\n");
  printf("      srfiYb = %7.1f +",srfiYorg);
  printf(" %7.4f * (srfiXb ",slope);
  if (srfiXorg>=0) then begin
     printf("- %6d)\n\n",srfiXorg);
  end
  else begin
     printf("+ %6d)\n\n",-srfiXorg);
  end
  printf("   SELECTED bnp VALUE: %5.2f\n\n",bnp);
  if (pVI == 1) then begin
     printf("   DENSE-VEGETATION POINT:\n");
     printf("      srfiXveg = %6d\n",srfiXv);
     printf("      srfiYveg = %6d\n\n",srfiYv);
  end
  else begin
     printf("   FOREGROUND-MATERIALS POINT:\n");
     printf("      srfiXfore = %6d\n",srfiXv);
     printf("      srfiYfore = %6d\n\n",srfiYv);
  end
  cvi = 1;
  if (pVI == 1) then begin
     if (tsrfi == 0) then begin
        if (abs(bnp) < 0.02) then begin
           printf("   RESULT: Classic NDVI\n");
           cvi = 0;
        end
        if (abs(bnp - 0.5) < 0.02) then begin
           printf("   RESULT: Classic SAVI\n");
           cvi = 0;
        end
        if (abs(bnp - 20) < 0.02) then begin
           printf("   RESULT: Classic PVI\n");
           cvi = 0;
        end
     end
     if (tsrfi == 1) then begin
        if (abs(bnp) < 0.02) then begin
           printf("   RESULT: Transformed NDVI\n");
           cvi = 0;
        end
        if (abs(bnp - 0.15) < 0.02) then begin
           printf("   RESULT: Classic TSAVI\n");
           cvi = 0;
        end
        if (abs(bnp - 20) < 0.02) then begin
           printf("   RESULT: Classic WDVI\n");
           cvi = 0;
        end
     end
     if (cvi == 1) then begin
        printf("   RESULT: Customized GRUVI\n");
     end
  end
  else begin
     printf("   RESULT: Customized GRUFI\n");
  end
  if (bnp < 0.01) then bnp = 0.01;
  bnpxy = bnp * 5000;
# ------------------------------------------------------------
# PERFORM TRANSLATION (T) & ROTATION (R) TRANSFORMATION
# TO CALCULATE GRUVI SCALE-CHANGE FACTOR, gfac:
  rang = -atan(slope);
  rangdeg = rang * deg; rangdeg = rangdeg + 45;
  sinrang = sind(rangdeg); cosrang = cosd(rangdeg);
  tsrfiX1 = srfiXv - srfiXorg;
  tsrfiY1 = srfiYv - srfiYorg;
  tsrfiX2 = tsrfiX1 * cosrang - tsrfiY1 * sinrang;
  tsrfiY2 = tsrfiX1 * sinrang + tsrfiY1 * cosrang;
  tsrfiX2 = tsrfiX2 + bnpxy;
  tsrfiY2 = tsrfiY2 + bnpxy;
  gruvidveg = (tsrfiY2 - tsrfiX2) / (tsrfiY2 + tsrfiX2); 
  gfac = gruvirange / gruvidveg;
  grubif = (bnpxy + srfiXorg) * (bnpxy + srfiXorg);
  grubif = grubif + (bnpxy + srfiYorg) * (bnpxy + srfiYorg);
  grubifmin = sqrt(grubif);
# ------------------------------------------------------------
# ASK USER TO GET INPUT RASTERS CALLED SRFIX & SRFIY:
  GetInputRaster(SRFIX);
  nlins = NumLins(SRFIX); ncols = NumCols(SRFIX);
  rtype$ = RastType(SRFIX);
  checkHisto(SRFIX);
  GetInputRaster(SRFIY,nlins,ncols,rtype$);
  checkHisto(SRFIY);
# ------------------------------------------------------------
# ASK USER TO DESIGNATE OUTPUT GRUVI & GRUBI RASTERS:
  GetOutputRaster(GRUBI,nlins,ncols,gtype$);
  printf("\nPROCESSING ACTIONS:\n");
  printf("Assigning GRUBI Values to Null Value.\n");
  GRUBI = grunull;
  SetNull(GRUBI,grunull);
  CopySubobjects(SRFIX,GRUBI,"GEOREF");
  if (pVI == 1) then begin
     GetOutputRaster(GRUVI,nlins,ncols,gtype$);
     printf("Assigning GRUVI Values to Null Value.\n");
     GRUVI = grunull;
     SetNull(GRUVI,grunull);
     CopySubobjects(SRFIX,GRUVI,"GEOREF");
  end
  else begin
     GetOutputRaster(GRUFI,nlins,ncols,gtype$);
     printf("Assigning GRUFI Values to Null Value.\n");
     GRUFI = grunull;
     SetNull(GRUFI,grunull);
     CopySubobjects(SRFIX,GRUFI,"GEOREF");
  end
  printf("Producing Output Rasters.\n");
  for each SRFIX begin
     isnull = IsNull(SRFIX);
     if (isnull == 0) then begin
        srfiX = SRFIX;
        srfiY = SRFIY;
        tsrfiX1 = srfiX - srfiXorg;
        tsrfiY1 = srfiY - srfiYorg;
        tsrfiX2 = tsrfiX1 * cosrang - tsrfiY1 * sinrang;
        tsrfiY2 = tsrfiX1 * sinrang + tsrfiY1 * cosrang;
        tsrfiX2 = tsrfiX2 + bnpxy;
        tsrfiY2 = tsrfiY2 + bnpxy;
        grubif = tsrfiX2 * tsrfiX2;
        grubif = grubif + tsrfiY2 * tsrfiY2;
        grubif = sqrt(grubif) - grubifmin;
        grubi = round(0.2 * grubif);
        gruvif = tsrfiY2 - tsrfiX2;
        gruvif = gfac * gruvif / (tsrfiY2 + tsrfiX2);
        gruvi = round(gruvif);
        if (grubi < -grumax) then grubi = -grumax;
        if (grubi > grumax) then grubi = grumax;
        if (gruvi < -grumax) then gruvi = -grumax;
        if (gruvi > grumax) then gruvi = grumax;
        GRUBI = grubi;
        if (pVI == 1) then begin
           GRUVI = gruvi;
        end
        else begin
           GRUFI = gruvi;
        end
     end
  end
# ------------------------------------------------------------
# CREATE HISTOGRAMS & PYRAMIDS.  THEN, CLOSE RASTERS:
  printf("Adding Histogram & Pyramid Subobjects.\n");
  if (pVI == 1) then begin
     CreateHistogram(GRUVI,0); CreatePyramid(GRUVI,0);
     hnmGF = HasNullMask(GRUVI); hnGF = HasNull(GRUVI);
     if (hnGF == 1) then nvGF = NullValue(GRUVI);
     c1GF = GlobalPtile(GRUVI,0.5,hnGF,nvGF);
     c3GF = GlobalPtile(GRUVI,50,hnGF,nvGF);
     c2GF = GlobalPtile(GRUVI,99.5,hnGF,nvGF);
     pGF = expPower(c1GF,c3GF,c2GF);
     smlContrast.InputLowerLimit = c1GF;
     smlContrast.InputUpperLimit = c2GF;
     smlContrast.Power           =  pGF;
     smlContrast.Compute(GRUVI,clut1$,clut2$,clut3$);
     CloseRaster(GRUVI);
  end
  else begin
     CreateHistogram(GRUFI,0); CreatePyramid(GRUFI,0);
     hnmGF = HasNullMask(GRUFI); hnGF = HasNull(GRUFI);
     if (hnGF == 1) then nvGF = NullValue(GRUFI);
     c1GF = GlobalPtile(GRUFI,0.5,hnGF,nvGF);
     c3GF = GlobalPtile(GRUFI,50,hnGF,nvGF);
     c2GF = GlobalPtile(GRUFI,99.5,hnGF,nvGF);
     pGF = expPower(c1GF,c3GF,c2GF);
     smlContrast.InputLowerLimit = c1GF;
     smlContrast.InputUpperLimit = c2GF;
     smlContrast.Power           =  pGF;
     smlContrast.Compute(GRUFI,clut1$,clut2$,clut3$);
     CloseRaster(GRUFI);
  end
  CreateHistogram(GRUBI,0); CreatePyramid(GRUBI);
  hnmGB = HasNullMask(GRUBI); hnGB = HasNull(GRUBI);
  if (hnGB == 1) then nvGB = NullValue(GRUBI);
  c1GB = GlobalPtile(GRUBI,0.5,hnGB,nvGB);
  c3GB = GlobalPtile(GRUBI,50,hnGB,nvGB);
  c2GB = GlobalPtile(GRUBI,99.5,hnGB,nvGB);
  pGB = expPower(c1GB,c3GB,c2GB);
  smlContrast.InputLowerLimit = c1GB;
  smlContrast.InputUpperLimit = c2GB;
  smlContrast.Power           =  pGB;
  smlContrast.Compute(GRUBI,clut1$,clut2$,clut3$);
  CloseRaster(SRFIX); CloseRaster(SRFIY);
  CloseRaster(GRUBI);
  printf("Processing Finished.");
  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.");