# ------------------------------------------------------------
# TASCAP.sml
# ------------------------------------------------------------
# SET WARNING LEVEL:
$warnings 3
# ------------------------------------------------------------
# DEFINE PROCEDURE writeTitle:
# PURPOSE: WRITES TITLE & AUTHOR INFO TO CONSOLE WINDOW.
proc writeTitle() begin
printf("TASCAP.sml:\n");
printf(" VERSION: November 15, 2005\n");
printf(" PURPOSE: CONVERT SRFI DATA ");
printf("TO TASSELED-CAP DATA\n");
printf(" WITH AUTOMATIC CONTRAST ");
printf("LOOK-UP TABLES.\n");
printf(" DETAILS: FAQs_by_Jack A & F\n");
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
# ------------------------------------------------------------
# DECLARE VARIABLES
numeric nbands,bp0,bp1,bp2,bp3,bp4,bp5,bp6,bp7,bp8,bp9;
numeric i,j,nbandsm1,m;
numeric itc,nTC,tc0,tc1,tc2,tc3,tc4,tc5,tc6,tc7,tc8,tc9;
numeric lin,col,nlins,ncols,fsize;
numeric tcn,tcnp1,isnull,lintcp,coltcp;
array lintcp[10]; array coltcp[10];
string bp0$,bp1$,bp2$,bp3$,bp4$,bp5$,bp6$,bp7$,bp8$;
string stype$,dtype$,ttype$;
dtype$ = "16-bit unsigned";
ttype$ = "16-bit signed"; tcn = -32000; tcnp1 = tcn + 1;
numeric d,d2,sd2,mag,f,sf2,ds,tc;
numeric dot1,dot2,dot3,dot4,dot5,dot6,dot7,dot8,dot9;
numeric c1DS0,c1DS1,c1DS2,c1DS3,c1DS4;
numeric c1DS5,c1DS6,c1DS7,c1DS8;
numeric c2DS0,c2DS1,c2DS2,c2DS3,c2DS4;
numeric c2DS5,c2DS6,c2DS7,c2DS8;
numeric c3DS0,c3DS1,c3DS2,c3DS3,c3DS4;
numeric c3DS5,c3DS6,c3DS7,c3DS8;
numeric pDS0,pDS1,pDS2,pDS3,pDS4,pDS5,pDS6,pDS7,pDS8;
numeric c1TC1,c1TC2,c1TC3,c1TC4,c1TC5,c1TC6,c1TC7,c1TC8;
numeric c2TC1,c2TC2,c2TC3,c2TC4,c2TC5,c2TC6,c2TC7,c2TC8;
numeric c3TC1,c3TC2,c3TC3,c3TC4,c3TC5,c3TC6,c3TC7,c3TC8;
numeric pTC1,pTC2,pTC3,pTC4,pTC5,pTC6,pTC7,pTC8;
numeric hnmDS0,hnmDS1,hnmDS2,hnmDS3,hnmDS4;
numeric hnmDS5,hnmDS6,hnmDS7,hnmDS8;
numeric hnDS0,hnDS1,hnDS2,hnDS3,hnDS4;
numeric hnDS5,hnDS6,hnDS7,hnDS8;
numeric nvDS0,nvDS1,nvDS2,nvDS3,nvDS4;
numeric nvDS5,nvDS6,nvDS7,nvDS8;
numeric hnmTC1,hnmTC2,hnmTC3,hnmTC4;
numeric hnmTC5,hnmTC6,hnmTC7,hnmTC8;
numeric hnTC1,hnTC2,hnTC3,hnTC4;
numeric hnTC5,hnTC6,hnTC7,hnTC8;
numeric nvTC1,nvTC2,nvTC3,nvTC4;
numeric nvTC5,nvTC6,nvTC7,nvTC8;
# ------------------------------------------------------------
# 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;
# ------------------------------------------------------------
# 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$;
# ------------------------------------------------------------
# DECLARE ARRAYS
array bp0[10]; array bp1[10]; array bp2[10]; array bp3[10];
array bp4[10]; array bp5[10]; array bp6[10]; array bp7[10];
array bp8[10]; array bp9[10]; array tc0[10]; array tc1[10];
array tc2[10]; array tc3[10]; array tc4[10]; array tc5[10];
array tc6[10]; array tc7[10]; array tc8[10]; array tc9[10];
array numeric srfi[10];
# ------------------------------------------------------------
# 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 METHOD (OUT OF 3 CHOICES): m
p1$ = "METHOD-NUMBER ENTRY:\n";
p2$ = " CHOICES:\n";
p3$ = " 1: Process LANDSAT SRFItoa DATA (6 Bands)\n";
p4$ = " Using USGS TASCAP Coefficients.\n";
p5$ = " 2: Process with PRE-DEFINED bp VALUES.\n";
p6$ = " 3: Specify KEY PIXELS (lin,col) to Fetch\n";
p7$ = " bp Values from Input SRFT Rasters.\n";
p8$ = "METHOD-NUMBER ENTERED:";
p$ = p1$ + p2$ + p3$ + p4$ + p5$ + p6$ + p7$ + p8$;
m = PopupNum(p$,1,1,3,0);
# ------------------------------------------------------------
# ASSIGN VALUES TO dtype$, ttype$, tcn, tcnp1, nbands & nTC:
dtype$ = "16-bit unsigned";
ttype$ = "16-bit signed"; tcn = -32000; tcnp1 = tcn + 1;
nbands = 6; nTC = 4;
# ============================================================
# DEFAULT PARAMETERS FOR METHOD 1: Landsat 6-Band MS Data.
# Classic TASCAP for SRFI Rasters: BL, GL, RL, NA, MB, MC:
# WARNING: Only Applicable to Landsat SRFItoa Raster Values.
if (m == 1) then begin
bp0$ = "BLACK"; bp1$ = "Brightness";
bp2$ = "Greenness" ; bp3$ = "Wetness" ;
bp4$ = "Haze" ;
tc0[1]= 0.0; tc0[2]= 0.0; tc0[3]= 0.0;
tc0[4]= 0.0; tc0[5]= 0.0; tc0[6]= 0.0;
tc1[1]= 0.3561; tc1[2]= 0.3972; tc1[3]= 0.3904;
tc1[4]= 0.6966; tc1[5]= 0.2286; tc1[6]= 0.1596;
tc2[1]= -0.3344; tc2[2]= -0.3544; tc2[3]= -0.4556;
tc2[4]= 0.6966; tc2[5]= -0.0242; tc2[6]= -0.2630;
tc3[1]= 0.2626; tc3[2]= 0.2141; tc3[3]= 0.0926;
tc3[4]= 0.0656; tc3[5]= -0.7629; tc3[6]= -0.5388;
tc4[1]= 0.0805; tc4[2]= -0.0498; tc4[3]= 0.1950;
tc4[4]= -0.1327; tc4[5]= 0.5752; tc4[6]= -0.7775;
end
# ============================================================
# DEFAULT PARAMETERS FOR METHOD 2:
# Default TASCAP for SRFI Rasters: BL, GL, RL, NA, (MB, MC):
# METHOD 2 Requires No Input Key Inputs from User:
# NOTE: METHOD 2 can be Applied to 3 to 6 Bands (SRFI):
if (m == 2) then begin
bp0$="Dark Soil "; bp1$="Bright Soil";
bp2$="Dense Green Veg."; bp3$="Water" ;
bp4$="Yellow Veg." ;
t$ = "n-SPACE ORIGIN TYPE";
p1$ = "SELECT n-SPACE ORIGIN TYPE:\n";
p2$ = " ENTER 'BLACK' or ACCEPT DEFAULT:\n";
p3$ = "n-SPACE ORIGIN-TYPE ENTERED:";
p$ = p1$ + p2$ + p3$;
bp0$ = PopupString(p$,bp0$,t$);
bp0[1]= 190; bp0[2]= 326; bp0[3]= 100;
bp0[4]= 200; bp0[5]=142; bp0[6]= 276;
if (bp0$ == "BLACK") then begin
bp0[1]=0; bp0[2]=0; bp0[3]=0;
bp0[4]=0; bp0[5]=0; bp0[6]=0;
end
bp1[1]= 945; bp1[2]=1157; bp1[3]=1746;
bp1[4]=1998; bp1[5]=2688; bp1[6]=1967;
bp2[1]= 400; bp2[2]= 600; bp2[3]= 300;
bp2[4]=6000; bp2[5]=2000; bp2[6]= 600;
bp3[1]= 400; bp3[2]= 500; bp3[3]= 200;
bp3[4]= 189; bp3[5]= 105; bp3[6]= 70;
bp4[1]= 450; bp4[2]=1000; bp4[3]=1400;
bp4[4]=4500; bp4[5]=2500; bp4[6]= 800;
end
# ============================================================
# METHOD 3: User must Provide Raster lin & col Locations that
# TASCAP.sml Will Use to Fetch Related bp Values.
# ============================================================
# 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
# ------------------------------------------------------------
# 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;
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 begin
rl = Rl;
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
# ------------------------------------------------------------
# WRITE TITLE & AUTHOR INFORMATION:
writeTitle();
# ------------------------------------------------------------
# WRITE METHOD TO CONSOLE WINDOW:
if (m == 1) then begin
printf("METHOD 1: Use USGS Coefficients\n");
printf(" Applied to Landsat SRFItoa Data.\n\n");
end
if (m == 2) then begin
printf("METHOD 2: Use Default TASCAP Coefficients\n");
printf(" Applied to SRFIsfc Data.\n\n");
end
if (m == 3) then begin
printf("METHOD 3: Use Customized TASCAP Coefficients\n");
printf(" Applied to SRFIsfc.\n");
printf(" Customized Names with Associated\n");
printf(" SRFI Image LIN & COL Locations.\n\n");
end
# ------------------------------------------------------------
# SPECIFY nbands & nTC PARAMETERS:
if (m == 2) then begin
p1$ = "NUMBER OF INPUT SRFI RASTERS ENTRY:\n";
p2$ = " RANGE: 4 to 6\n";
p3$ = "NUMBER ENTERED:";
p$ = p1$ + p2$ + p3$;
nbands = PopupNum(p$,nbands,4,6,0);
nbandsm1 = nbands - 1;
end
if (m == 2) then begin
if (nTC > (nbands - 1)) then nTC = nbands - 1;
p1$ = "NUMBER OF OUTPUT RASTER PAIRS ENTRY:\n";
p2$ = " RANGE: 2 to (NumInputBands - 1)\n";
p3$ = "NUMBER ENTERED:";
p$ = p1$ + p2$ + p3$;
nTC = PopupNum(p$,nTC,2,nbandsm1,0);
end
if (m == 3) then begin
p1$ = "NUMBER OF INPUT SRFI RASTERS ENTRY:\n";
p2$ = " RANGE: 3 to 9\n";
p3$ = "NUMBER ENTERED:";
p$ = p1$ + p2$ + p3$;
nbands = PopupNum(p$,nbands,3,9,0);
nbandsm1 = nbands - 1;
end
printf("NUMBERS OF INPUT RASTERS & OUTPUT RASTERS:\n");
printf(" %1d Sets of SRFI INPUT Rasters\n",nbands);
if (m == 3) then begin
if (nTC > (nbands - 1)) then nTC = nbands - 1;
p1$ = "NUMBER OF OUTPUT RASTER PAIRS ENTRY:\n";
p2$ = " RANGE: 2 to (NumInputBands - 1)\n";
p3$ = "NUMBER ENTERED:";
p$ = p1$ + p2$ + p3$;
nTC = PopupNum(p$,nTC,2,nbandsm1,0);
end
printf(" %1d Pairs of TC & DS OUTPUT Rasters\n\n",nTC);
# ------------------------------------------------------------
# OPEN INPUT RASTERS:
raster SRFI1,SRFI2,SRFI3,SRFI4,SRFI5;
raster SRFI6,SRFI7,SRFI8,SRFI9;
GetInputRaster(SRFI1);
nlins = NumLins(SRFI1); ncols = NumCols(SRFI1);
stype$ = RastType(SRFI1);
checkHisto(SRFI1);
GetInputRaster(SRFI2,nlins,ncols,stype$);
checkHisto(SRFI2);
GetInputRaster(SRFI3,nlins,ncols,stype$);
checkHisto(SRFI3);
if (nbands>3) then begin
GetInputRaster(SRFI4,nlins,ncols,stype$);
checkHisto(SRFI4);
end
if (nbands>4) then begin
GetInputRaster(SRFI5,nlins,ncols,stype$);
checkHisto(SRFI5);
end
if (nbands>5) then begin
GetInputRaster(SRFI6,nlins,ncols,stype$);
checkHisto(SRFI6);
end
if (nbands>6) then begin
GetInputRaster(SRFI7,nlins,ncols,stype$);
checkHisto(SRFI7);
end
if (nbands>7) then begin
GetInputRaster(SRFI8,nlins,ncols,stype$);
checkHisto(SRFI8);
end
if (nbands>8) then begin
GetInputRaster(SRFI9,nlins,ncols,stype$);
checkHisto(SRFI9);
end
# ------------------------------------------------------------
# GET USER-PROVIDED LIN & COL PARAMETERS:
if (m == 3) then begin
t$ = "METHOD-3 PARAMETER";
bp0$ = "bp0 NAME";
p$ = "ENTER bp0 NAME:";
bp0$ = PopupString(p$,bp0$,t$);
if (bp0$ == "BLACK") then begin
lintcp[0] = 1;
coltcp[0] = 1;
end
p1$ = "METHOD-3 LINE & COLUMN PARAMETERS ENTRY:\n";
if (bp0$ <> "BLACK") then begin
p2$ = "FOR " + bp0$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[0] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp0$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[0] = PopupNum(p$,1,1,100000,0);
end
bp1$ = "bp1 NAME";
p$ = "ENTER bp1 NAME:";
bp1$ = PopupString(p$,bp1$,t$);
p2$ = "FOR " + bp1$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[1] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp1$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[1] = PopupNum(p$,1,1,100000,0);
bp2$ = "bp2 NAME";
p$ = "ENTER bp2 NAME:";
bp2$ = PopupString(p$,bp2$,t$);
p2$ = "FOR " + bp2$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[2] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp2$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[2] = PopupNum(p$,1,1,100000,0);
if (nTC > 2) then begin
bp3$ = "bp3 NAME";
p$ = "ENTER bp3 NAME:";
bp3$ = PopupString(p$,bp3$,t$);
p2$ = "FOR " + bp3$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[3] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp3$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[3] = PopupNum(p$,1,1,100000,0);
end
if (nTC > 3) then begin
bp4$ = "bp4 NAME";
p$ = "ENTER bp4 NAME:";
bp4$ = PopupString(p$,bp4$,t$);
p2$ = "FOR " + bp4$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[4] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp4$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[4] = PopupNum(p$,1,1,100000,0);
end
if (nTC > 4) then begin
bp5$ = "bp5 NAME";
p$ = "ENTER bp5 NAME:";
bp5$ = PopupString(p$,bp5$,t$);
p2$ = "FOR " + bp5$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[5] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp5$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[5] = PopupNum(p$,1,1,100000,0);
end
if (nTC > 5) then begin
bp6$ = "bp6 NAME";
p$ = "ENTER bp6 NAME:";
bp6$ = PopupString(p$,bp6$,t$);
p2$ = "FOR " + bp6$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[6] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp6$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[6] = PopupNum(p$,1,1,100000,0);
end
if (nTC > 6) then begin
bp7$ = "bp7 NAME";
p$ = "ENTER bp7 NAME:";
bp7$ = PopupString(p$,bp7$,t$);
p2$ = "FOR " + bp7$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[7] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp7$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[7] = PopupNum(p$,1,1,100000,0);
end
if (nTC > 7) then begin
bp8$ = "bp8 NAME";
p$ = "ENTER bp8 NAME:";
bp8$ = PopupString(p$,bp8$,t$);
p2$ = "FOR " + bp8$ + ", ENTER LIN:";
p$ = p1$ + p2$;
lintcp[8] = PopupNum(p$,1,1,100000,0);
p2$ = "FOR " + bp8$ + ", ENTER COL:";
p$ = p1$ + p2$;
coltcp[8] = PopupNum(p$,1,1,100000,0);
end
end
# ------------------------------------------------------------
# GET bpn[nband] VALUES:
for i=0 to nTC begin
if (m == 3) then begin
lin = lintcp[i]; col = coltcp[i];
if (i==0) then begin
bp0[1] = 0; bp0[2] = 0; bp0[3] = 0;
bp0[4] = 0; bp0[5] = 0; bp0[6] = 0;
bp0[7] = 0; bp0[8] = 0; bp0[9] = 0;
if (bp0$ <> "BLACK") then begin
bp0[1] = SRFI1[lin,col]; bp0[2] = SRFI2[lin,col];
bp0[3] = SRFI3[lin,col];
if (nbands > 3) then bp0[4] = SRFI4[lin,col];
if (nbands > 4) then bp0[5] = SRFI5[lin,col];
if (nbands > 5) then bp0[6] = SRFI6[lin,col];
if (nbands > 6) then bp0[7] = SRFI7[lin,col];
if (nbands > 7) then bp0[8] = SRFI8[lin,col];
if (nbands > 8) then bp0[9] = SRFI9[lin,col];
end
end
if (i==1) then begin
bp1[1] = SRFI1[lin,col]; bp1[2] = SRFI2[lin,col];
bp1[3] = SRFI3[lin,col];
if (nbands > 3) then bp1[4] = SRFI4[lin,col];
if (nbands > 4) then bp1[5] = SRFI5[lin,col];
if (nbands > 5) then bp1[6] = SRFI6[lin,col];
if (nbands > 6) then bp1[7] = SRFI7[lin,col];
if (nbands > 7) then bp1[8] = SRFI8[lin,col];
if (nbands > 8) then bp1[9] = SRFI9[lin,col];
end
if (i==2) then begin
bp2[1] = SRFI1[lin,col]; bp2[2] = SRFI2[lin,col];
bp2[3] = SRFI3[lin,col];
if (nbands > 3) then bp2[4] = SRFI4[lin,col];
if (nbands > 4) then bp2[5] = SRFI5[lin,col];
if (nbands > 5) then bp2[6] = SRFI6[lin,col];
if (nbands > 6) then bp2[7] = SRFI7[lin,col];
if (nbands > 7) then bp2[8] = SRFI8[lin,col];
if (nbands > 8) then bp2[9] = SRFI9[lin,col];
end
if (i==3) then begin
bp3[1] = SRFI1[lin,col]; bp3[2] = SRFI2[lin,col];
bp3[3] = SRFI3[lin,col];
if (nbands > 3) then bp3[4] = SRFI4[lin,col];
if (nbands > 4) then bp3[5] = SRFI5[lin,col];
if (nbands > 5) then bp3[6] = SRFI6[lin,col];
if (nbands > 6) then bp3[7] = SRFI7[lin,col];
if (nbands > 7) then bp3[8] = SRFI8[lin,col];
if (nbands > 8) then bp3[9] = SRFI9[lin,col];
end
if (i==4) then begin
bp4[1] = SRFI1[lin,col]; bp4[2] = SRFI2[lin,col];
bp4[3] = SRFI3[lin,col];
if (nbands > 3) then bp4[4] = SRFI4[lin,col];
if (nbands > 4) then bp4[5] = SRFI5[lin,col];
if (nbands > 5) then bp4[6] = SRFI6[lin,col];
if (nbands > 6) then bp4[7] = SRFI7[lin,col];
if (nbands > 7) then bp4[8] = SRFI8[lin,col];
if (nbands > 8) then bp4[9] = SRFI9[lin,col];
end
if (i==5) then begin
bp5[1] = SRFI1[lin,col]; bp5[2] = SRFI2[lin,col];
bp5[3] = SRFI3[lin,col];
if (nbands > 3) then bp5[4] = SRFI4[lin,col];
if (nbands > 4) then bp5[5] = SRFI5[lin,col];
if (nbands > 5) then bp5[6] = SRFI6[lin,col];
if (nbands > 6) then bp5[7] = SRFI7[lin,col];
if (nbands > 7) then bp5[8] = SRFI8[lin,col];
if (nbands > 8) then bp5[9] = SRFI9[lin,col];
end
if (i==6) then begin
bp6[1] = SRFI1[lin,col]; bp6[2] = SRFI2[lin,col];
bp6[3] = SRFI3[lin,col];
if (nbands > 3) then bp6[4] = SRFI4[lin,col];
if (nbands > 4) then bp6[5] = SRFI5[lin,col];
if (nbands > 5) then bp6[6] = SRFI6[lin,col];
if (nbands > 6) then bp6[7] = SRFI7[lin,col];
if (nbands > 7) then bp6[8] = SRFI8[lin,col];
if (nbands > 8) then bp6[9] = SRFI9[lin,col];
end
if (i==7) then begin
bp7[1] = SRFI1[lin,col]; bp7[2] = SRFI2[lin,col];
bp7[3] = SRFI3[lin,col];
if (nbands > 3) then bp7[4] = SRFI4[lin,col];
if (nbands > 4) then bp7[5] = SRFI5[lin,col];
if (nbands > 5) then bp7[6] = SRFI6[lin,col];
if (nbands > 6) then bp7[7] = SRFI7[lin,col];
if (nbands > 7) then bp7[8] = SRFI8[lin,col];
if (nbands > 8) then bp7[9] = SRFI9[lin,col];
end
if (i==8) then begin
bp8[1] = SRFI1[lin,col]; bp8[2] = SRFI2[lin,col];
bp8[3] = SRFI3[lin,col];
if (nbands > 3) then bp8[4] = SRFI4[lin,col];
if (nbands > 4) then bp8[5] = SRFI5[lin,col];
if (nbands > 5) then bp8[6] = SRFI6[lin,col];
if (nbands > 6) then bp8[7] = SRFI7[lin,col];
if (nbands > 7) then bp8[8] = SRFI8[lin,col];
if (nbands > 8) then bp8[9] = SRFI9[lin,col];
end
if (i==9) then begin
bp9[1] = SRFI1[lin,col]; bp9[2] = SRFI2[lin,col];
bp9[3] = SRFI3[lin,col];
if (nbands > 3) then bp9[4] = SRFI4[lin,col];
if (nbands > 4) then bp9[5] = SRFI5[lin,col];
if (nbands > 5) then bp9[6] = SRFI6[lin,col];
if (nbands > 6) then bp9[7] = SRFI7[lin,col];
if (nbands > 7) then bp9[8] = SRFI8[lin,col];
if (nbands > 8) then bp9[9] = SRFI9[lin,col];
end
end
end
# ------------------------------------------------------------
# WRTIE BIOPHYSICAL VECTORS TO CONSOLE WINDOW (SRFI Units):
printf("BIOPHYSICAL CONTROL MATERIALS:\n");
printf(" BPO: %s\n",bp0$);
if (m > 1) then begin
printf(" BP0:");
nbandsm1 = nbands - 1;
for i=1 to nbandsm1 begin
printf("%8d",bp0[i]);
end
printf("%8d\n\n",bp0[nbands]);
end
printf(" BP1: %s\n",bp1$);
if (m > 1) then begin
printf(" BP1:");
for i=1 to nbandsm1 begin
printf("%8d",bp1[i]);
end
printf("%8d\n\n",bp1[nbands]);
end
printf(" BP2: %s\n",bp2$);
if (m > 1) then begin
printf(" BP2:");
for i=1 to nbandsm1 begin
printf("%8d",bp2[i]);
end
printf("%8d\n\n",bp2[nbands]);
end
for itc=3 to nTC begin
if (itc==3) then begin
printf(" BP3: %s\n",bp3$);
if (m > 1) then begin
printf(" BP3:");
for i=1 to nbandsm1 begin
printf("%8d",bp3[i]);
end
printf("%8d\n\n",bp3[nbands]);
end
end
if (itc==4) then begin
printf(" BP4: %s\n",bp4$);
if (m > 1) then begin
printf(" BP4:");
for i=1 to nbandsm1 begin
printf("%8d",bp4[i]);
end
printf("%8d\n\n",bp4[nbands]);
end
end
if (itc==5) then begin
printf(" BP5: %s\n",bp5$);
if (m > 1) then begin
printf(" BP5:");
for i=1 to nbandsm1 begin
printf("%8d",bp5[i]);
end
printf("%8d\n\n",bp5[nbands]);
end
end
if (itc==6) then begin
printf(" BP6: %s\n",bp6$);
if (m > 1) then begin
printf(" BP6:");
for i=1 to nbandsm1 begin
printf("%8d",bp6[i]);
end
printf("%8d\n\n",bp6[nbands]);
end
end
if (itc==7) then begin
printf(" BP7: %s\n",bp7$);
if (m > 1) then begin
printf(" BP7:");
for i=1 to nbandsm1 begin
printf("%8d",bp7[i]);
end
printf("%8d\n\n",bp7[nbands]);
end
end
if (itc==8) then begin
printf(" BP8: %s\n",bp8$);
if (m > 1) then begin
printf(" BP8:");
for i=1 to nbandsm1 begin
printf("%8d",bp8[i]);
end
printf("%8d\n\n",bp8[nbands]);
end
end
end
# ------------------------------------------------------------
# COMPUTE TC OFFSETS & ROTATION COEFFICIENTS:
printf("TC OFFSETS (tc0) COEFFICIENTS (tc1 etc.):\n");
# tc0:
printf(" tc0: ");
for i=1 to nbands begin
if (m > 1) then begin
tc0[i] = bp0[i];
end
printf("%8d",tc0[i]);
end
# tc1:
printf("\n tc1: ");
if (m > 1) then begin
sd2 = 0;
for i=1 to nbands begin
d = bp1[i] - tc0[i];
sd2 = sd2 + d * d;
end
mag = sqrt(sd2);
end
for i=1 to nbands begin
if (m > 1) then begin
d = bp1[i] - tc0[i];
tc1[i] = d / mag;
end
printf("%8.4f",tc1[i]);
end
# tc2:
printf("\n tc2: ");
if (m > 1) then begin
dot1 = 0;
for i=1 to nbands begin
d = bp2[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp2[i] - tc0[i];
f = d - dot1 * tc1[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
end
for i=1 to nbands begin
if (m > 1) then begin
d = bp2[i] - tc0[i];
f = d - dot1 * tc1[i];
tc2[i] = f / mag;
end
printf("%8.4f",tc2[i]);
end
# tc3:
if (nTC > 2) then begin
printf("\n tc3: ");
if (m > 1) then begin
dot1 = 0; dot2 = 0;
for i=1 to nbands begin
d = bp3[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
dot2 = dot2 + tc2[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp3[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
end
for i=1 to nbands begin
if (m > 1) then begin
d = bp3[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
tc3[i] = f / mag;
end
printf("%8.4f",tc3[i]);
end
end
# tc4:
if (nTC > 3) then begin
printf("\n tc4: ");
if (m > 1) then begin
dot1 = 0; dot2 = 0; dot3 = 0;
for i=1 to nbands begin
d = bp4[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
dot2 = dot2 + tc2[i] * d;
dot3 = dot3 + tc3[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp4[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
end
for i=1 to nbands begin
if (m > 1) then begin
d = bp4[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
tc4[i] = f / mag;
end
printf("%8.4f",tc4[i]);
end
end
# tc5:
if (nTC > 4) then begin
printf("\n tc5: ");
dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0;
for i=1 to nbands begin
d = bp5[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
dot2 = dot2 + tc2[i] * d;
dot3 = dot3 + tc3[i] * d;
dot4 = dot4 + tc4[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp5[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
for i=1 to nbands begin
d = bp5[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
tc5[i] = f / mag;
printf("%8.4f",tc5[i]);
end
end
# tc6:
if (nTC > 5) then begin
printf("\n tc6: ");
dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
for i=1 to nbands begin
d = bp6[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
dot2 = dot2 + tc2[i] * d;
dot3 = dot3 + tc3[i] * d;
dot4 = dot4 + tc4[i] * d;
dot5 = dot5 + tc5[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp6[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
for i=1 to nbands begin
d = bp6[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
tc6[i] = f / mag;
printf("%8.4f",tc6[i]);
end
end
# tc7:
if (nTC > 6) then begin
printf("\n tc7: ");
dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
dot6 = 0;
for i=1 to nbands begin
d = bp7[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
dot2 = dot2 + tc2[i] * d;
dot3 = dot3 + tc3[i] * d;
dot4 = dot4 + tc4[i] * d;
dot5 = dot5 + tc5[i] * d;
dot6 = dot6 + tc6[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp7[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
f = f - dot6 * tc6[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
for i=1 to nbands begin
d = bp7[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
f = f - dot6 * tc6[i];
tc7[i] = f / mag;
printf("%8.4f",tc7[i]);
end
end
# tc8:
if (nTC > 7) then begin
printf("\n tc8: ");
dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
dot6 = 0; dot7 = 0;
for i=1 to nbands begin
d = bp8[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
dot2 = dot2 + tc2[i] * d;
dot3 = dot3 + tc3[i] * d;
dot4 = dot4 + tc4[i] * d;
dot5 = dot5 + tc5[i] * d;
dot6 = dot6 + tc6[i] * d;
dot7 = dot7 + tc7[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp8[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
f = f - dot6 * tc6[i];
f = f - dot7 * tc7[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
for i=1 to nbands begin
d = bp8[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
f = f - dot6 * tc6[i];
f = f - dot7 * tc7[i];
tc8[i] = f / mag;
printf("%8.4f",tc8[i]);
end
end
# tc9:
if (nTC > 6) then begin
printf("\n tc9: ");
dot1 = 0; dot2 = 0; dot3 = 0; dot4 = 0; dot5 = 0;
dot6 = 0; dot7 = 0; dot8 = 0;
for i=1 to nbands begin
d = bp9[i] - tc0[i];
dot1 = dot1 + tc1[i] * d;
dot2 = dot2 + tc2[i] * d;
dot3 = dot3 + tc3[i] * d;
dot4 = dot4 + tc4[i] * d;
dot5 = dot5 + tc5[i] * d;
dot6 = dot6 + tc6[i] * d;
dot7 = dot7 + tc7[i] * d;
dot8 = dot8 + tc8[i] * d;
end
sf2 = 0;
for i=1 to nbands begin
d = bp9[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
f = f - dot6 * tc6[i];
f = f - dot7 * tc7[i];
f = f - dot8 * tc8[i];
sf2 = sf2 + f * f;
end
mag = sqrt(sf2);
for i=1 to nbands begin
d = bp9[i] - tc0[i];
f = d - dot1 * tc1[i];
f = f - dot2 * tc2[i];
f = f - dot3 * tc3[i];
f = f - dot4 * tc4[i];
f = f - dot5 * tc5[i];
f = f - dot6 * tc6[i];
f = f - dot7 * tc7[i];
f = f - dot8 * tc8[i];
tc9[i] = f / mag;
printf("%8.4f",tc9[i]);
end
end
# ------------------------------------------------------------
# SETUP OUTPUT RASTERS:
raster DS0,TC1,DS1,TC2,DS2,TC3,DS3,TC4,DS4,TC5,DS5;
raster TC6,DS6,TC7,DS7,TC8,DS8;
GetOutputRaster(DS0,nlins,ncols,dtype$); SetNull(DS0,0);
CopySubobjects(SRFI1,DS0,"GEOREF");
GetOutputRaster(TC1,nlins,ncols,ttype$); SetNull(TC1,tcn);
CopySubobjects(SRFI1,TC1,"GEOREF");
GetOutputRaster(DS1,nlins,ncols,dtype$); SetNull(DS1,0);
CopySubobjects(SRFI1,DS1,"GEOREF");
GetOutputRaster(TC2,nlins,ncols,ttype$); SetNull(TC2,tcn);
CopySubobjects(SRFI1,TC2,"GEOREF");
GetOutputRaster(DS2,nlins,ncols,dtype$); SetNull(DS2,0);
CopySubobjects(SRFI1,DS2,"GEOREF");
if (nTC>2) then begin
GetOutputRaster(TC3,nlins,ncols,ttype$);
SetNull(TC3,tcn);
CopySubobjects(SRFI1,TC3,"GEOREF");
GetOutputRaster(DS3,nlins,ncols,dtype$); SetNull(DS3,0);
CopySubobjects(SRFI1,DS3,"GEOREF");
end
if (nTC>3) then begin
GetOutputRaster(TC4,nlins,ncols,ttype$);
SetNull(TC4,tcn);
CopySubobjects(SRFI1,TC4,"GEOREF");
GetOutputRaster(DS4,nlins,ncols,dtype$); SetNull(DS4,0);
CopySubobjects(SRFI1,DS4,"GEOREF");
end
if (nTC>4) then begin
GetOutputRaster(TC5,nlins,ncols,ttype$);
SetNull(TC5,tcn);
CopySubobjects(SRFI1,TC5,"GEOREF");
GetOutputRaster(DS5,nlins,ncols,dtype$); SetNull(DS5,0);
CopySubobjects(SRFI1,DS5,"GEOREF");
end
if (nTC>5) then begin
GetOutputRaster(TC6,nlins,ncols,ttype$);
SetNull(TC6,tcn);
CopySubobjects(SRFI1,TC6,"GEOREF");
GetOutputRaster(DS6,nlins,ncols,dtype$); SetNull(DS6,0);
CopySubobjects(SRFI1,DS6,"GEOREF");
end
if (nTC>6) then begin
GetOutputRaster(TC7,nlins,ncols,ttype$);
SetNull(TC7,tcn);
CopySubobjects(SRFI1,TC7,"GEOREF");
GetOutputRaster(DS7,nlins,ncols,dtype$); SetNull(DS7,0);
CopySubobjects(SRFI1,DS7,"GEOREF");
end
if (nTC>7) then begin
GetOutputRaster(TC8,nlins,ncols,ttype$);
SetNull(TC8,tcn);
CopySubobjects(SRFI1,TC8,"GEOREF");
GetOutputRaster(DS8,nlins,ncols,dtype$); SetNull(DS8,0);
CopySubobjects(SRFI1,DS8,"GEOREF");
end
# ------------------------------------------------------------
# COMPUTE DS & TC VALUES & ASSIGN TO DS & TC RASTERS:
printf("\n\nPRODUCING OUTPUT RASTERS:\n");
printf(" DS0");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
sd2 = 0;
for i=1 to nbands begin
d = srfi[i] - tc0[i];
sd2 = sd2 + d * d;
end
ds = round(sqrt(sd2));
if (ds < 1) then ds = 1;
DS0 = ds;
end
end
printf(" TC1 DS1");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS0;
tc = 0;
for i=1 to nbands begin
tc = tc + tc1[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC1 = tc; DS1 = ds;
end
end
printf(" TC2 DS2");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS1; tc = 0;
for i=1 to nbands begin
tc = tc + tc2[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC2 = tc; DS2 = ds;
end
end
if (nTC>2) then begin
printf(" TC3 DS3");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS2; tc = 0;
for i=1 to nbands begin
tc = tc + tc3[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC3 = tc; DS3 = ds;
end
end
end
if (nTC>3) then begin
printf(" TC4 DS4");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS3; tc = 0;
for i=1 to nbands begin
tc = tc + tc4[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC4 = tc; DS4 = ds;
end
end
end
if (nTC>4) then begin
printf(" TC5 DS5");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS4; tc = 0;
for i=1 to nbands begin
tc = tc + tc5[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC5 = tc; DS5 = ds;
end
end
end
if (nTC>5) then begin
printf(" TC6 DS6");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS5; tc = 0;
for i=1 to nbands begin
tc = tc + tc6[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC6 = tc; DS6 = ds;
end
end
end
if (nTC>6) then begin
printf(" TC7 DS7");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS6; tc = 0;
for i=1 to nbands begin
tc = tc + tc7[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC7 = tc; DS7 = ds;
end
end
end
if (nTC>7) then begin
printf(" TC8 DS8");
for each SRFI1 begin
isnull = IsNull(SRFI1);
if (isnull <> 1) then begin
srfi[1] = SRFI1; srfi[2] = SRFI2; srfi[3] = SRFI3;
if (nbands>3) then srfi[4] = SRFI4;
if (nbands>4) then srfi[5] = SRFI5;
if (nbands>5) then srfi[6] = SRFI6;
if (nbands>6) then srfi[7] = SRFI7;
if (nbands>7) then srfi[8] = SRFI8;
if (nbands>8) then srfi[9] = SRFI9;
ds = DS7; tc = 0;
for i=1 to nbands begin
tc = tc + tc8[i] * (srfi[i] - tc0[i]);
end
ds = ds * ds - tc * tc;
if (ds < 1) then ds = 1;
ds = round(sqrt(ds));
if (tc < tcnp1) then tc = tcnp1;
TC8 = tc; DS8 = ds;
end
end
end
printf("\n\n");
# ------------------------------------------------------------
# CLOSE INPUT RASTERS:
CloseRaster(SRFI1); CloseRaster(SRFI2); CloseRaster(SRFI3);
if (nbands>3) then CloseRaster(SRFI4);
if (nbands>4) then CloseRaster(SRFI5);
if (nbands>5) then CloseRaster(SRFI6);
if (nbands>6) then CloseRaster(SRFI7);
if (nbands>7) then CloseRaster(SRFI8);
if (nbands>8) then CloseRaster(SRFI9);
# ------------------------------------------------------------
# MAKE SUBOBJECTS:
printf("MAKING SUBOBJECTS:\n");
printf("CONTRAST LOOK-UP TABLE PARAMETERS:\n");
printf("RASTER LOW HIGH POWER\n");
CreateHistogram(DS0,0); CreatePyramid(DS0);
hnmDS0 = HasNullMask(DS0); hnDS0 = HasNull(DS0);
if (hnDS0 == 1) then nvDS0 = NullValue(DS0);
c1DS0 = GlobalPtile(DS0,0.5,hnDS0,nvDS0);
c3DS0 = GlobalPtile(DS0,50,hnDS0,nvDS0);
c2DS0 = GlobalPtile(DS0,99.5,hnDS0,nvDS0);
pDS0 = expPower(c1DS0,c3DS0,c2DS0);
smlContrast.InputLowerLimit = c1DS0;
smlContrast.InputUpperLimit = c2DS0;
smlContrast.Power = pDS0;
smlContrast.Compute(DS0,clut1$,clut2$,clut3$);
printf(" DS0 %9d %9d %9.2f\n",c1DS0,c2DS0,pDS0);
CloseRaster(DS0);
CreateHistogram(TC1,0); CreatePyramid(TC1);
hnmTC1 = HasNullMask(TC1); hnTC1 = HasNull(TC1);
if (hnTC1 == 1) then nvTC1 = NullValue(TC1);
c1TC1 = GlobalPtile(TC1,0.5,hnTC1,nvTC1);
c3TC1 = GlobalPtile(TC1,50,hnTC1,nvTC1);
c2TC1 = GlobalPtile(TC1,99.5,hnTC1,nvTC1);
pTC1 = expPower(c1TC1,c3TC1,c2TC1);
smlContrast.InputLowerLimit = c1TC1;
smlContrast.InputUpperLimit = c2TC1;
smlContrast.Power = pTC1;
smlContrast.Compute(TC1,clut1$,clut2$,clut3$);
printf(" TC1 %9d %9d %9.2f\n",c1TC1,c2TC1,pTC1);
CloseRaster(TC1);
CreateHistogram(DS1,0); CreatePyramid(DS1);
hnmDS1 = HasNullMask(DS1); hnDS1 = HasNull(DS1);
if (hnDS1 == 1) then nvDS1 = NullValue(DS1);
c1DS1 = GlobalPtile(DS1,0.5,hnDS1,nvDS1);
c3DS1 = GlobalPtile(DS1,50,hnDS1,nvDS1);
c2DS1 = GlobalPtile(DS1,99.5,hnDS1,nvDS1);
pDS1 = expPower(c1DS1,c3DS1,c2DS1);
smlContrast.InputLowerLimit = c1DS1;
smlContrast.InputUpperLimit = c2DS1;
smlContrast.Power = pDS1;
smlContrast.Compute(DS1,clut1$,clut2$,clut3$);
printf(" DS1 %9d %9d %9.2f\n",c1DS1,c2DS1,pDS1);
CloseRaster(DS1);
CreateHistogram(TC2,0); CreatePyramid(TC2);
hnmTC2 = HasNullMask(TC2); hnTC2 = HasNull(TC2);
if (hnTC2 == 1) then nvTC2 = NullValue(TC2);
c1TC2 = GlobalPtile(TC2,0.5,hnTC2,nvTC2);
c3TC2 = GlobalPtile(TC2,50,hnTC2,nvTC2);
c2TC2 = GlobalPtile(TC2,99.5,hnTC2,nvTC2);
pTC2 = expPower(c1TC2,c3TC2,c2TC2);
smlContrast.InputLowerLimit = c1TC2;
smlContrast.InputUpperLimit = c2TC2;
smlContrast.Power = pTC2;
smlContrast.Compute(TC2,clut1$,clut2$,clut3$);
printf(" TC2 %9d %9d %9.2f\n",c1TC2,c2TC2,pTC2);
CloseRaster(TC2);
CreateHistogram(DS2,0); CreatePyramid(DS2);
hnmDS2 = HasNullMask(DS2); hnDS2 = HasNull(DS2);
if (hnDS2 == 1) then nvDS2 = NullValue(DS2);
c1DS2 = GlobalPtile(DS2,0.5,hnDS2,nvDS2);
c3DS2 = GlobalPtile(DS2,50,hnDS2,nvDS2);
c2DS2 = GlobalPtile(DS2,99.5,hnDS2,nvDS2);
pDS2 = expPower(c1DS2,c3DS2,c2DS2);
smlContrast.InputLowerLimit = c1DS2;
smlContrast.InputUpperLimit = c2DS2;
smlContrast.Power = pDS2;
smlContrast.Compute(DS2,clut1$,clut2$,clut3$);
printf(" DS2 %9d %9d %9.2f\n",c1DS2,c2DS2,pDS2);
CloseRaster(DS2);
if (nTC>2) then begin
CreateHistogram(TC3,0); CreatePyramid(TC3);
hnmTC3 = HasNullMask(TC3); hnTC3 = HasNull(TC3);
if (hnTC3 == 1) then nvTC3 = NullValue(TC3);
c1TC3 = GlobalPtile(TC3,0.5,hnTC3,nvTC3);
c3TC3 = GlobalPtile(TC3,50,hnTC3,nvTC3);
c2TC3 = GlobalPtile(TC3,99.5,hnTC3,nvTC3);
pTC3 = expPower(c1TC3,c3TC3,c2TC3);
smlContrast.InputLowerLimit = c1TC3;
smlContrast.InputUpperLimit = c2TC3;
smlContrast.Power = pTC3;
smlContrast.Compute(TC3,clut1$,clut2$,clut3$);
printf(" TC3 %9d %9d %9.2f\n",c1TC3,c2TC3,pTC3);
CloseRaster(TC3);
CreateHistogram(DS3,0); CreatePyramid(DS3);
hnmDS3 = HasNullMask(DS3); hnDS3 = HasNull(DS3);
if (hnDS3 == 1) then nvDS3 = NullValue(DS3);
c1DS3 = GlobalPtile(DS3,0.5,hnDS3,nvDS3);
c3DS3 = GlobalPtile(DS3,50,hnDS3,nvDS3);
c2DS3 = GlobalPtile(DS3,99.5,hnDS3,nvDS3);
pDS3 = expPower(c1DS3,c3DS3,c2DS3);
smlContrast.InputLowerLimit = c1DS3;
smlContrast.InputUpperLimit = c2DS3;
smlContrast.Power = pDS3;
smlContrast.Compute(DS3,clut1$,clut2$,clut3$);
printf(" DS3 %9d %9d %9.2f\n",c1DS3,c2DS3,pDS3);
CloseRaster(DS3);
end
if (nTC>3) then begin
CreateHistogram(TC4,0); CreatePyramid(TC4);
hnmTC4 = HasNullMask(TC4); hnTC4 = HasNull(TC4);
if (hnTC4 == 1) then nvTC4 = NullValue(TC4);
c1TC4 = GlobalPtile(TC4,0.5,hnTC4,nvTC4);
c3TC4 = GlobalPtile(TC4,50,hnTC4,nvTC4);
c2TC4 = GlobalPtile(TC4,99.5,hnTC4,nvTC4);
pTC4 = expPower(c1TC4,c3TC4,c2TC4);
smlContrast.InputLowerLimit = c1TC4;
smlContrast.InputUpperLimit = c2TC4;
smlContrast.Power = pTC4;
smlContrast.Compute(TC4,clut1$,clut2$,clut3$);
printf(" TC4 %9d %9d %9.2f\n",c1TC4,c2TC4,pTC4);
CloseRaster(TC4);
CreateHistogram(DS4,0); CreatePyramid(DS4);
hnmDS4 = HasNullMask(DS4); hnDS4 = HasNull(DS4);
if (hnDS4 == 1) then nvDS4 = NullValue(DS4);
c1DS4 = GlobalPtile(DS4,0.5,hnDS4,nvDS4);
c3DS4 = GlobalPtile(DS4,50,hnDS4,nvDS4);
c2DS4 = GlobalPtile(DS4,99.5,hnDS4,nvDS4);
pDS4 = expPower(c1DS4,c3DS4,c2DS4);
smlContrast.InputLowerLimit = c1DS4;
smlContrast.InputUpperLimit = c2DS4;
smlContrast.Power = pDS4;
smlContrast.Compute(DS4,clut1$,clut2$,clut3$);
printf(" DS4 %9d %9d %9.2f\n",c1DS4,c2DS4,pDS4);
CloseRaster(DS4);
end
if (nTC>4) then begin
CreateHistogram(TC5,0); CreatePyramid(TC5);
hnmTC5 = HasNullMask(TC5); hnTC5 = HasNull(TC5);
if (hnTC5 == 1) then nvTC5 = NullValue(TC5);
c1TC5 = GlobalPtile(TC5,0.5,hnTC5,nvTC5);
c3TC5 = GlobalPtile(TC5,50,hnTC5,nvTC5);
c2TC5 = GlobalPtile(TC5,99.5,hnTC5,nvTC5);
pTC5 = expPower(c1TC5,c3TC5,c2TC5);
smlContrast.InputLowerLimit = c1TC5;
smlContrast.InputUpperLimit = c2TC5;
smlContrast.Power = pTC5;
smlContrast.Compute(TC5,clut1$,clut2$,clut3$);
printf(" TC5 %9d %9d %9.2f\n",c1TC5,c2TC5,pTC5);
CloseRaster(TC5);
CreateHistogram(DS5,0); CreatePyramid(DS5);
hnmDS5 = HasNullMask(DS5); hnDS5 = HasNull(DS5);
if (hnDS5 == 1) then nvDS5 = NullValue(DS5);
c1DS5 = GlobalPtile(DS5,0.5,hnDS5,nvDS5);
c3DS5 = GlobalPtile(DS5,50,hnDS5,nvDS5);
c2DS5 = GlobalPtile(DS5,99.5,hnDS5,nvDS5);
pDS5 = expPower(c1DS5,c3DS5,c2DS5);
smlContrast.InputLowerLimit = c1DS5;
smlContrast.InputUpperLimit = c2DS5;
smlContrast.Power = pDS5;
smlContrast.Compute(DS5,clut1$,clut2$,clut3$);
printf(" DS5 %9d %9d %9.2f\n",c1DS5,c2DS5,pDS5);
CloseRaster(DS5);
end
if (nTC>5) then begin
CreateHistogram(TC6,0); CreatePyramid(TC6);
hnmTC6 = HasNullMask(TC6); hnTC6 = HasNull(TC6);
if (hnTC6 == 1) then nvTC6 = NullValue(TC6);
c1TC6 = GlobalPtile(TC6,0.5,hnTC6,nvTC6);
c3TC6 = GlobalPtile(TC6,50,hnTC6,nvTC6);
c2TC6 = GlobalPtile(TC6,99.5,hnTC6,nvTC6);
pTC6 = expPower(c1TC6,c3TC6,c2TC6);
smlContrast.InputLowerLimit = c1TC6;
smlContrast.InputUpperLimit = c2TC6;
smlContrast.Power = pTC6;
smlContrast.Compute(TC6,clut1$,clut2$,clut3$);
printf(" TC6 %9d %9d %9.2f\n",c1TC6,c2TC6,pTC6);
CloseRaster(TC6);
CreateHistogram(DS6,0); CreatePyramid(DS6);
hnmDS6 = HasNullMask(DS6); hnDS6 = HasNull(DS6);
if (hnDS6 == 1) then nvDS6 = NullValue(DS6);
c1DS6 = GlobalPtile(DS6,0.5,hnDS6,nvDS6);
c3DS6 = GlobalPtile(DS6,50,hnDS6,nvDS6);
c2DS6 = GlobalPtile(DS6,99.5,hnDS6,nvDS6);
pDS6 = expPower(c1DS6,c3DS6,c2DS6);
smlContrast.InputLowerLimit = c1DS6;
smlContrast.InputUpperLimit = c2DS6;
smlContrast.Power = pDS6;
smlContrast.Compute(DS6,clut1$,clut2$,clut3$);
printf(" DS6 %9d %9d %9.2f\n",c1DS6,c2DS6,pDS6);
CloseRaster(DS6);
end
if (nTC>6) then begin
CreateHistogram(TC7,0); CreatePyramid(TC7);
hnmTC7 = HasNullMask(TC7); hnTC7 = HasNull(TC7);
if (hnTC7 == 1) then nvTC7 = NullValue(TC7);
c1TC7 = GlobalPtile(TC7,0.5,hnTC7,nvTC7);
c3TC7 = GlobalPtile(TC7,50,hnTC7,nvTC7);
c2TC7 = GlobalPtile(TC7,99.5,hnTC7,nvTC7);
pTC7 = expPower(c1TC7,c3TC7,c2TC7);
smlContrast.InputLowerLimit = c1TC7;
smlContrast.InputUpperLimit = c2TC7;
smlContrast.Power = pTC7;
smlContrast.Compute(TC7,clut1$,clut2$,clut3$);
printf(" TC7 %9d %9d %9.2f\n",c1TC7,c2TC7,pTC7);
CloseRaster(TC7);
CreateHistogram(DS7,0); CreatePyramid(DS7);
hnmDS7 = HasNullMask(DS7); hnDS7 = HasNull(DS7);
if (hnDS7 == 1) then nvDS7 = NullValue(DS7);
c1DS7 = GlobalPtile(DS7,0.5,hnDS7,nvDS7);
c3DS7 = GlobalPtile(DS7,50,hnDS7,nvDS7);
c2DS7 = GlobalPtile(DS7,99.5,hnDS7,nvDS7);
pDS7 = expPower(c1DS7,c3DS7,c2DS7);
smlContrast.InputLowerLimit = c1DS7;
smlContrast.InputUpperLimit = c2DS7;
smlContrast.Power = pDS7;
smlContrast.Compute(DS7,clut1$,clut2$,clut3$);
printf(" DS7 %9d %9d %9.2f\n",c1DS7,c2DS7,pDS7);
CloseRaster(DS7);
end
if (nTC>7) then begin
CreateHistogram(TC8,0); CreatePyramid(TC8);
hnmTC8 = HasNullMask(TC8); hnTC8 = HasNull(TC8);
if (hnTC8 == 1) then nvTC8 = NullValue(TC8);
c1TC8 = GlobalPtile(TC8,0.5,hnTC8,nvTC8);
c3TC8 = GlobalPtile(TC8,50,hnTC8,nvTC8);
c2TC8 = GlobalPtile(TC8,99.5,hnTC8,nvTC8);
pTC8 = expPower(c1TC8,c3TC8,c2TC8);
smlContrast.InputLowerLimit = c1TC8;
smlContrast.InputUpperLimit = c2TC8;
smlContrast.Power = pTC8;
smlContrast.Compute(TC8,clut1$,clut2$,clut3$);
printf(" TC8 %9d %9d %9.2f\n",c1TC8,c2TC8,pTC8);
CloseRaster(TC8);
CreateHistogram(DS8,0); CreatePyramid(DS8);
hnmDS8 = HasNullMask(DS8); hnDS8 = HasNull(DS8);
if (hnDS8 == 1) then nvDS8 = NullValue(DS8);
c1DS8 = GlobalPtile(DS8,0.5,hnDS8,nvDS8);
c3DS8 = GlobalPtile(DS8,50,hnDS8,nvDS8);
c2DS8 = GlobalPtile(DS8,99.5,hnDS8,nvDS8);
pDS8 = expPower(c1DS8,c3DS8,c2DS8);
smlContrast.InputLowerLimit = c1DS8;
smlContrast.InputUpperLimit = c2DS8;
smlContrast.Power = pDS8;
smlContrast.Compute(DS8,clut1$,clut2$,clut3$);
printf(" DS8 %9d %9d %9.2f\n",c1DS8,c2DS8,pDS8);
CloseRaster(DS8);
end
printf("\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.");
Download