# ------------------------------------------------------------
# WATER.sml
# ------------------------------------------------------------
# SET WARNING LEVEL:
  $warnings 3
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle:
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
  proc writeTitle() begin
     printf("WATER.sml:\n");
     printf("         VERSION: November 15, 2005\n");
     printf("         PURPOSE: MAKE A MERGED WATER & ");
     printf("LAND IMAGE.\n");
     printf("         DETAILS: FAQs_by_Jack G\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
# ------------------------------------------------------------
# DECLARE VARIABLES:
  numeric nlins,ncols;
  numeric m,d,tLP,v,nc,wm,pW,pL,x,y,yt,facBL,offNA;
  numeric sBL_8,sGL_8,sRL_8,sNA_8;
  numeric srfiBL,srfiNA;
  array sBL_8[65536],sGL_8[65536],sRL_8[65536],sNA_8[65536];
  numeric i,mini,srfiBLmaxW,srfiGLmaxW,srfiRLmaxW;
  numeric srfiBLmaxL,srfiGLmaxL,srfiRLmaxL,srfiNAmaxL;
  numeric redderW,brighterW,brightenDarkW;
  numeric redderL,brighterL,brightenDarkL;
  numeric pW,pL;
  string stype$;
  raster SRFIBL,SRFIGL,SRFIRL,SRFINA,WATERMASK;
  raster IMAGE;
  numeric iB,iG,iR,fB,fG,fR,fCLR;
# ------------------------------------------------------------
# DEFINE VALUES FOR fR, and fG:
  fR = 2^16; fG = 2^8; fB = 1;
# ------------------------------------------------------------
# DEFINE STRING VARIABLES FOR POPUP WINDOWS:
  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$;
# ------------------------------------------------------------
# CLEAR CONSOLE WINDOW & REQUEST REPOSITIONING:
  p1$ = "CONSOLE-WINDOW ADJUSTMENT:\n";
  p2$ = "* REPOSITION the CONSOLE WINDOW.\n";
  p3$ = "* Then, CLICK OK.";
  p$  = p1$ + p2$ + p3$;
  PopupMessage(p$);
# ------------------------------------------------------------
# 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
# ------------------------------------------------------------
# USER INPUTS:
# ENHANCEDMENT DEFAULTS:
  redderW       = 1.18;
  brighterW     = 1.00;
  brightenDarkW = 1.43;
  redderL       = 0.95;
  brighterL     = 1.00;
  brightenDarkL = 1.25;
# ------------------------------------------------------------
# TNTmips VERSION / PATCH INPUT:
  p1$ = "VERSION / PATCH ENTRY:\n";
  p2$ = "   Yes: Running TNTmips V7.1 With PATCH\n";
  p3$ = "        Dated July 13, 2005, or LATER\n";
  p4$ = "    No: Running Earlier VERSION or PATCH.";
  p$ = p1$ + p2$ + p3$ + p4$;
  fCLR = PopupYesNo(p$,1);
  if (fCLR == 1) then begin
     fR = 1; fB = 2^16;
  end
# ------------------------------------------------------------
# ACCEPT ENHANCEMENT DEFAULTS:
  p$ = "Do You Want to Apply Default ENHANCEMENTS?";
  d =  PopupYesNo(p$,0);
# ------------------------------------------------------------
# WATER-MASK OPTION:
  p1$ = "WATER-MASK OPTION ENTRY:\n";
  p2$ = "   OPTIONS:\n";
  p3$ = "      1: Create a NEW WATERMASK Raster.\n";
  p4$ = "      2: Use an EXISTING WATERMASK Raster.\n";
  p5$ = "         Existing Raster Must Match SRFI Rasters.\n";
  p6$ = "OPTION ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ + p6$;
  m = PopupNum(p$,2,1,2,0);
# ------------------------------------------------------------
# LAND-IMAGE OPTION:
  p1$ = "LAND-IMAGE OPTION ENTRY:\n";
  p2$ = "   OPTIONS:\n";
  p3$ = "      1: COLOR INFRARED IMAGE (CIR).\n";
  p4$ = "      2. NATURAL COLOR IMAGE.\n";
  p5$ = "      3: GRAYSCALE (GREEN BAND ONLY) IMAGE.\n";
  p6$ = "      4: GRAYSCALE (NIR BAND ONLY) IMAGE.\n";
  p7$ = "OPTION ENTERED:";
  p$  = p1$ + p2$ + p3$ + p4$ + p5$ +p6$ +p7$;
  tLP = PopupNum(p$,4,1,4,0);
# ------------------------------------------------------------
# MAKE WATER-COLORS REDDER:
  p1$ = "REDDER-WATER-COLORS PARAMETER ENTRY:\n"; 
  p2$ = "VALUE OF redderW ENTERED:";
  p$ = p1$ + p2$;
  if (d == 0) then begin
     redderW = PopupNum(p$,redderW,0.1,2,2);
  end
# ------------------------------------------------------------
# MAKE ALL-WATER-COLORS BRIGHTER:
  p1$ = "BRIGHTER-WATER-COLORS PARAMETER ENTRY:\n";
  p2$ = "VALUE FOR brighterW ENTERED:";
  p$ = p1$ + p2$;
  if (d == 0) then begin
     brighterW = PopupNum(p$,brighterW,0.1,2,2);
  end
# ------------------------------------------------------------
# MAKE DARK-WATER-COLORS BRIGHTER:
  p1$ = "BRIGHTEN-DARK-WATER-COLORS PARAMETER ENTRY:\n";
  p2$ = "VALUE OF brightenDarkW ENTERED:";
  p$ = p1$ + p2$;
  if (d == 0) then begin
     brightenDarkW = PopupNum(p$,brightenDarkW,0.1,2,2);
  end
# ------------------------------------------------------------
# MAKE LAND-COLORS REDDER:
  p1$ = "REDDER-LAND-COLORS PARAMETER ENTRY:\n";
  p2$ = "VALUE OF redderL ENTERED:";
  p$ = p1$ + p2$;
  if (d == 0) then begin
     if (tLP < 3) then begin
        redderL = PopupNum(p$,redderL,0.1,2,2);
     end
  end
# ------------------------------------------------------------
# MAKE ALL-LAND-COLORS BRIGHTER:
  p1$ = "BRIGHTER-LAND-COLORS PARAMETER ENTRY:\n";
  p2$ = "VALUE OF brighterL ENTERED:";
  p$ = p1$ + p2$;
  if (d == 0) then begin
     brighterL = PopupNum(p$,brighterL,0.1,2,2);
  end
# ------------------------------------------------------------
# MAKE DARK-LAND-COLORS BRIGHTER:
  p1$ = "BRIGHTEN-DARK-LAND-COLORS PARAMETER ENTRY:\n";
  p2$ = "VALUE OF brightenDarkL ENTERED:";
  p$ = p1$ + p2$;
  if (d == 0) then begin
     brightenDarkL = PopupNum(p$,brightenDarkL,0.1,2,2);
  end
# ------------------------------------------------------------
# LAND-WATER BOUNDARY CONTROL PARAMETERS:
  if (m == 1) then begin
     p1$ = "LINE-EQUATION MULTIPLIER ENTRY:\n";
     p2$ = "   RANGE: 0.00 to 10.00\n";
     p3$ = "VALUE OF facBL ENTERED:";
     p$ = p1$ + p2$ + p3$;
     facBL = PopupNum(p$,0.16,0.0,10,2);
     p1$ = "LINE-EQUATION SRFINA OFFSET ENTRY:\n";
     p2$ = "   RANGE: 1 to 2000\n";
     p3$ = "VALUE OF offNA ENTERED:";
     p$ = p1$ + p2$ + p3$;
     offNA = PopupNum(p$,925,1,3000,0);
  end
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION:
  clear();
  writeTitle();
# ------------------------------------------------------------
# PRINT USER INPUTS TO CONSOLE WINDOW:
  if (m == 1) then begin
     printf("A NEW WATER-MASK IS BEING CREATED.\n\n");
  end
  if (m == 2) then begin
     printf("AN EXISTING WATER-MASK IS BEING USED.\n");
     printf("NOTE: Existing WATER-MASK Must Match");
     printf(" SRFI rasters.\n\n");
  end
  if (tLP == 1) then begin
     printf("LAND PIXELS WILL HAVE CIR COLORS.");
  end
  if (tLP == 2) then begin
     printf("LAND PIXELS WILL HAVE NATURAL COLORS.");
  end
  if (tLP == 3) then begin
     printf("LAND PIXELS WILL HAVE GREEN-BAND GRAY TONES.");
  end
  if (tLP == 4) then begin
     printf("LAND PIXELS WILL HAVE NIR-BAND GRAY TONES.");
  end
  printf("\n\n");
  printf("IMAGE ENHANCEMENT PARAMETERS FOR WATER PIXELS:\n");
  printf("         redderW = %4.2f\n",redderW);
  printf("       brighterW = %4.2f\n",brighterW);
  printf("   brightenDarkW = %4.2f\n\n",brightenDarkW);
  printf("IMAGE ENHANCEMENT PARAMETERS FOR LAND PIXELS:\n");
  if (tLP < 3) then begin
     printf("         redderL = %4.2f\n",redderL);
  end
  printf("       brighterL = %4.2f\n",brighterL);
  printf("   brightenDarkL = %4.2f\n",brightenDarkL);
  printf("\n");
  if (m == 1) then begin
     printf("WATERMASK CONTROL PARAMETERS:\n");
     printf("        facBL = %4.2f\n",facBL);
     printf("        offNA = %4d\n\n",offNA);
  end
# ------------------------------------------------------------
# OPEN SRFIsfc RASTERS:
  printf("OPEN SRFIsfc RASTERS (ANY APPROPRIATE SOURCE):\n");
  printf("   WARNING!  MUST be SURFACE-RELATED.\n");
  printf("   SRFIBL");
  GetInputRaster(SRFIBL);
  nlins = NumLins(SRFIBL); ncols = NumCols(SRFIBL);
  stype$ = RastType(SRFIBL);
  checkHisto(SRFIBL);
  printf("  SRFIGL");
  GetInputRaster(SRFIGL,nlins,ncols,stype$);
  checkHisto(SRFIGL);
  printf("  SRFIRL");
  GetInputRaster(SRFIRL,nlins,ncols,stype$);
  checkHisto(SRFIRL);
  printf("  SRFINA");
  GetInputRaster(SRFINA,nlins,ncols,stype$);
  checkHisto(SRFINA);
  printf("\n\n");
# ------------------------------------------------------------
# IF m == 2, OPEN WATER-MASK RASTER:
  if (m == 2) then begin
     printf("OPEN WATERMASK RASTER.");
     GetInputRaster(WATERMASK,nlins,ncols,"binary");
     checkHisto(WATERMASK);
     IgnoreNull(WATERMASK);
     printf("\n\n");
  end
# ------------------------------------------------------------
# IF m == 1, CREATE NEW OUTPUT RASTER(S:
  printf("CREATE NEW OUTPUT RASTER(S):\n");
  if (m == 1) then begin
     printf("   WATERMASK");
     GetOutputRaster(WATERMASK,nlins,ncols,"binary");
     IgnoreNull(WATERMASK);
     CopySubobjects(SRFIBL,WATERMASK,"GEOREF");
     printf("  IMAGE");
     GetOutputRaster(IMAGE,nlins,ncols,"24-bit color RGB");
     IgnoreNull(IMAGE);
     CopySubobjects(SRFIBL,IMAGE,"GEOREF");
  end
# ------------------------------------------------------------
# IF m == 2, CREATE NEW IMAGE RASTER ONLY:
  if (m == 2) then begin
     printf("   IMAGE");
     GetOutputRaster(IMAGE,nlins,ncols,"24-bit color RGB");
     IgnoreNull(IMAGE);
     CopySubobjects(SRFIBL,IMAGE,"GEOREF");
  end
  printf("\n\n");
# ------------------------------------------------------------
# PRODUCE VALUES FOR OUTPUT WATERMASK PIXELS:
  if (m == 1) then begin
     printf("FILLING WATERMASK RASTER WITH VALUES.");
        for each SRFINA begin
           x = SRFIBL; nc = IsNull(x);
           if (nc == 0) then begin
              wm = 1; yt = round(facBL * x + offNA);
              y = SRFINA;
              if (y > yt) then wm = 0;
           end
           WATERMASK = wm;
        end
     printf("\n\n");
  end
# ------------------------------------------------------------
# PRODUCE VALUES FOR OUTPUT IMAGE RASTER WATER-PIXELS:
  srfiGLmaxW = 1500 / brighterW;
  srfiRLmaxW = srfiGLmaxW / (redderW^1.51);
  srfiBLmaxW = srfiGLmaxW * redderW;
  pW = 1 / brightenDarkW;
  printf("DERIVED ENHANCEMENT FACTORS FOR WATER PIXELS:\n");
  printf("   srfiBLmaxW = %4d\n",srfiBLmaxW);
  printf("   srfiGLmaxW = %4d\n",srfiGLmaxW);
  printf("   srfiRLmaxW = %4d\n",srfiRLmaxW);
  printf("           pW = %4.2f\n\n",pW);
  for i=1 to 65535 begin
     v = 1 + round(254*(((i - 1) / srfiBLmaxW)^pW));
     if (v < 1) then v = 1;
     if (v > 255) then v = 255;
     sBL_8[i] = v;
  end
  for i=1 to 65535 begin
     v = 1 + round(254*(((i - 1) / srfiGLmaxW)^pW));
     if (v < 1) then v = 1;
     if (v > 255) then v = 255;
     sGL_8[i] = v;
  end
  for i=1 to 65535 begin
     v = 1 + round(254*(((i - 1) / srfiRLmaxW)^pW));
     if (v < 1) then v = 1;
     if (v > 255) then v = 255;
     sRL_8[i] = v;
  end
  if (tLP == 1) then begin
     srfiRLmaxL = 5500 / brighterL;
     srfiNAmaxL = srfiRLmaxL / redderL;
     srfiGLmaxL = srfiRLmaxL * (redderL^1.01);
     srfiBLmaxL = srfiRLmaxL * (redderL^1.68);
  end
  if (tLP == 2) then begin
     srfiGLmaxL = 3951 / brighterL;
     srfiBLmaxL = srfiGLmaxL * redderL;
     srfiRLmaxL = srfiGLmaxL / (redderL^1.65);
     srfiNAmaxL = srfiGLmaxL / (redderL^1.97);
  end
  if (tLP == 3) then begin
     srfiGLmaxL = 3951 / brighterL;
     srfiBLmaxL = srfiGLmaxL * redderL;
     srfiRLmaxL = srfiGLmaxL / (redderL^1.65);
     srfiNAmaxL = srfiGLmaxL / (redderL^1.97);
  end
  if (tLP == 4) then begin
     srfiNAmaxL = 5800 / brighterL;
     srfiRLmaxL = srfiNAmaxL * redderL;
     srfiGLmaxL = srfiNAmaxL * (redderL^1.99);
     srfiBLmaxL = srfiNAmaxL * (redderL^2.64);
  end
  pL = 1 / brightenDarkL;
  printf("DERIVED ENHANCEMENT FACTORS FOR LAND PIXELS:\n");
  if (tLP == 2) then begin
     printf("   srfiBLmaxL = %4d\n",srfiBLmaxL);
  end
  if (tLP < 4) then begin
     printf("   srfiGLmaxL = %4d\n",srfiGLmaxL);
  end
  if (tLP < 3) then begin
     printf("   srfiRLmaxL = %4d\n",srfiRLmaxL);
  end
  if (tLP == 1 or tLP == 4) then begin
     printf("   srfiNAmaxL = %4d\n",srfiNAmaxL);
  end
  printf("           pL = %4.2f\n\n",pL);
  printf("COLORIZING THE WATER PIXELS.");
  for each WATERMASK begin
     wm = WATERMASK;
     if (wm == 1) then begin
        v = SRFIBL; iB = sBL_8[v]; v = SRFIGL;
        iG = sGL_8[v]; v = SRFIRL; iR= sRL_8[v];
        IMAGE = iR * fR + iG * fG + iB * fB; 
     end
  end
  printf("\n\n");
# ------------------------------------------------------------
# PRODUCE VALUES FOR OUTPUT IMAGE RASTERS LAND-PIXELS:
  for i=1 to 65535 begin
     v = 1 + round(254*(((i - 1) / srfiBLmaxL)^pL));
     if (v < 1) then v = 1;
     if (v > 255) then v = 255;
     sBL_8[i] = v;
  end
  for i=1 to 65535 begin
     v = 20 + round(235*(((i - 1) / srfiGLmaxL)^pL));
     if (v < 1) then v = 1;
     if (v > 255) then v = 255;
     sGL_8[i] = v;
  end
  for i=1 to 65535 begin
     v = 1 + round(254*(((i - 1) / srfiRLmaxL)^pL));
     if (v < 1) then v = 1;
     if (v > 255) then v = 255;
     sRL_8[i] = v;
  end
  for i=1 to 65535 begin
     v = 1 + round(254*(((i - 1) / srfiNAmaxL)^pL));
     if (v < 1) then v = 1;
     if (v > 255) then v = 255;
     sNA_8[i] = v;
  end
  printf("COLORIZING THE LAND PIXELS.");
  for each WATERMASK begin
     wm = WATERMASK;
     if (wm == 0) then begin
        v = SRFIGL; nc = IsNull(v);
        if (nc == 0) then begin
           if (tLP <> 2) then begin
              iB = sGL_8[v]; v = SRFIRL; iG = sRL_8[v];
              v = SRFINA; iR = sNA_8[v];
           end
           if (tLP == 2) then begin
              v = SRFIBL; iB = sBL_8[v]; v = SRFIGL;
              iG = sGL_8[v]; v = SRFIRL; iR = sRL_8[v];
           end
           if (tLP == 3) then begin
              iG = iB; iR = iB;
           end
           if (tLP == 4) then begin
              iB = iR; iG = iR;
           end
           IMAGE = iR * fR + iG * fG + iB * fB;
        end
     end
  end
  printf("\n\n");
  printf("CREATING HISTOGRAMS & PYRAMIDS FOR IMAGE RASTERS.");
# ------------------------------------------------------------
# CREATE HISTOGRAMS & PYRAMIDS FOR OUTPUT RASTERS:
  SetNull(IMAGE,0); CreateHistogram(IMAGE,0);
  CreatePyramid(IMAGE,0);
  if (m == 1) then begin
     CreateHistogram(WATERMASK,0); CreatePyramid(WATERMASK,0);
  end
# ------------------------------------------------------------
# CLOSE or DELETE RASTERS:
  CloseRaster(SRFIBL); CloseRaster(SRFIGL);
  CloseRaster(SRFIRL); CloseRaster(SRFINA);
  CloseRaster(WATERMASK); CloseRaster(IMAGE);
  printf("\n\n");
  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.");