#### deveg68.sml
#### Sample script to subdue the expression of vegetation and enhance
#### the expression of underlying rock and soil in a multispectral
#### satellite image. Requires near-infrared and red bands to
#### determine the vegetation expression. Vegetation can be suppressed
#### in the near-infrared and red bands and in any additional
#### image bands selected. All input rasters must have the same
#### line-column dimensions and cell size.
#### This script implements a method developed by Robert E. Crippen
#### and Ronald G. Blom of NASA Jet Propulsion Laboratory:
#### Crippen, R.E. and Blom, R.G., 2001, Unveiling the Lithology of
#### Vegetated Terrains in Remotely Sensed Imagery, Photogrammetric
#### Engineering and Remote Sensing, Vol. 67, No. 8, pp. 935-943.
#### Products created by the script include:
#### a Vegetation Index raster ( VI )
#### and for each input band to be devegetated
#### a devegetated band raster
#### a 2DCorr raster showing a scatterplot of band values
#### versus VI
#### a Graph CAD object with a line representing the smoothed
#### curve fit to the plot of band value versus VI
#### Prior to running the script, determine a dark-pixel (path
#### radiance) value for each input raster to subtract effects of
#### sensor gain and brightness due to atmospheric scattering (haze).
#### Use the "Set dark pixel values..." button to enter these values.
#### (The default value is the minimum value for the band).
#### Dialog windows opened by the script are constructed using dialog
#### specifications in XML embedded in the script as multi-line string
#### variables.
#### This script requires TNTmips 6.8 release version or later.
#### Version of devegX.sml revised to deal with change in counter
#### behavior for class STRINGLIST adopted for TNTmips 6.8. (The previous
#### DevegX.sml no longer runs properly in TNTmips 6.8). Also corrected error
#### in retrieval of average band value for each VI level, which resulted in
#### output null values for VI = 0.
#### Randall Smith, MicroImages, Inc., 12 June 2003.
######################################
######## Global Variable Declarations
######################################
class XMLDOC docmain; # class instance for XML text containing main dialog specification
class XMLNODE nodemain; # XMLNODE for main dialog window
class GUI_DLG dlgmain; # dialog class for main dialog window
class XMLDOC docdp; # class instance for XML text contain dark-pixel correction dialog
class XMLNODE nodedp; # XMLNODE for dark-pixel correction dialog
class GUI_DLG dlgdp; # dialog class for dark-pixel correction dialog
string xmlmain$; # string variable for main dialog specification
numeric errXML, dlgreturn; # variables for error checking with main XML and dialog
class GUI_FORMDATA winsettings; # structure to hold all values set in main dialog window
class RVC_RASTER NIR, RED; # input near-infrared and red image bands
numeric lins, cols; # number of lines and columns in input raster set
string datatype$; # raster type for input image bands
string nirfile$, nirobjname$; # filename and object name for input NIR band
string redfile$, redobjname$; # filename and object name for input RED band
numeric nirobjnum, redobjnum; # object numbers for input NIR and RED bands
string nirtext$, redtext$; # strings with filename and object name for input NIR and RED
numeric devegNIR, devegRED; # flags to indicate if NIR and Red bands should
# be devegetated.
numeric numbands; # number of additional bands to be devegetated
numeric num; # counter for loops to process additional input bands
class RVC_RASTER Rin; # raster variable for additional input bands
string infile$, inobjname$; # name of file and object selected for Rin
numeric inobjnum; # object number of Rin in its parent file
numeric nirmin, redmin; # minimum values for input NIR and RED bands
# used as defaults for dark-pixel correction
numeric prNIR, prRED; # user-defined dark-pixel values to be subtracted
string prompt$; # prompt string for pop-up dialog to select dark-pixel-correction value
string infilepath$, infilename$; # components of input file name string
string deffile$; # default output Project File name from input file name
string outfile$; # name of output Project File
string mipspath$; # directory path to current TNTmips installation
string cmapfilename$; # path and filename for file containing standard color palette objects
numeric cmapnum; # object number of color palette object to be copied to correlation raster
class RVC_RASTER NIRlookup, REDlookup, PRlookup; # temporary rasters for lookup values
# for dark-pixel-correction
class RVC_RASTER RawVI; # temporary raster for raw vegetation index
class RVC_RASTER VI; # output raster for rescaled integer vegetation index
numeric n, r; # dark-pixel-corrected value read from lookup table for NIR, RED;
numeric minRawVI, maxRawVI; # min and max values for raw vegetation index
numeric range; # range of raw RawVI
numeric uplim, lowlim; # upper and lower RawVI limits to be used for rescaling
# to integer VI
numeric scale, offset; # scale and offset values used for rescaling RawVI to
# integer VI
class RVC_RASTER Corr2D, Corr2Dnir, Corr2Dred; # raster versions of 2D scatterplots of VI versus band value
CAD GraphNIR, GraphRED, Graph; # CAD objects for smoothed VI versus band value plots
numeric initavg, smoothcurve; # flags set by user to set smoothing parameters for
# fitting curve to the VI versus Band Value scatterplot
class RVC_RASTER NIRdeveg, REDdeveg; # rasters for devegetated NIR and RED bands
string nirdvobjname$, reddvobjname$; # object names for devegetated NIR and RED bands
numeric num; # counter for loop to devegetate additional bands
string num$; # string version of counter for status message
string corr$, corrdesc$; # strings for object name and description for correlation rasters
# for additional devegetated bands
class RVC_RASTER Deveg; # raster variable for output devegetated additional band
string outobjname$; # object name for output devegetated additional band
string cadname$; # object name for smoothed VI versus band value graph
#####################################################################
#####################################################################
################ Procedure definitions ############################
#####################################################################
#####################################################################
#####################################
### Callback procedure to initally disable OK button
### on main dialog window when it opens
#####################################
proc OnOpenMain () {
dlgmain.SetOkEnabled( 0 );
}
##################################
### Procedure to select NIR band
##################################
proc SelectNIR () {
GetInputRaster( NIR );
nirmin = GlobalMin( NIR );
# printf( "NIR min = %d\n", nirmin );
lins = NumLins( NIR );
cols = NumCols( NIR );
datatype$ = RastType( NIR );
nirfile$ = GetObjectFileName( NIR ); # directory path and input project file name
nirobjnum = GetObjectNumber( NIR );
nirobjname$ = GetObjectName( nirfile$, nirobjnum );
nirtext$ = nirfile$ + " \ " + nirobjname$;
dlgmain.SetCtrlValueStr( "NIRtext", nirtext$ );
dlgmain.GetCtrlByID( "getred" ).SetEnabled( 1 );
} # end SelectNIR ()
##################################
### Procedure to select RED band
##################################
proc SelectRED () {
GetInputRaster( RED, lins, cols, datatype$ );
redmin = GlobalMin( RED );
# printf( "RED min = %d\n", redmin );
redfile$ = GetObjectFileName( RED ); # directory path and input project file name
redobjnum = GetObjectNumber( RED );
redobjname$ = GetObjectName( redfile$, redobjnum );
redtext$ = redfile$ + " \ " + redobjname$;
dlgmain.SetCtrlValueStr("REDtext",redtext$);
dlgmain.GetCtrlByID("setdp").SetEnabled(1);
} # end SelectRED ()
##################################
### Callback procedure when number of bands
### is changed
###################################
proc OnNumChng () {
numbands = dlgmain.GetCtrlByID("numbands").GetValueNum();
if ( numbands > 0 ) {
dlgmain.GetCtrlByID("getaddbands").SetEnabled(1);
dlgmain.UpdateLayout();
}
}
#####################################
### Procedure to select additional bands for processing
#####################################
proc GetAddBands () {
if ( numbands > 0 ) {
class STRINGLIST inrastfilelist, inrastobjlist; # lists of strings with input raster
# filename(s) and object names
array numeric rastmins[ numbands ]; # array with minimum values for additional input rasters
local string listitem$;
class GUI_CTRL_LISTBOX addbandlist;
addbandlist = dlgmain.GetCtrlByID("bandlist");
addbandlist.SetEnabled(1);
for num = 1 to numbands {
GetInputRaster( Rin, lins, cols, datatype$ );
rastmins[ num ] = GlobalMin( Rin );
# printf("Min for add band %d = %d\n", num, rastmins[ num ] );
inobjnum = GetObjectNumber( Rin );
infile$ = GetObjectFileName( Rin );
inobjname$ = GetObjectName( infile$, inobjnum );
inrastfilelist.AddToEnd( infile$ );
inrastobjlist.AddToEnd( inobjname$ );
listitem$ = infile$ + " \ " + inobjname$;
addbandlist.AddItem( listitem$ );
}
}
} # end GetAddBands ()
#######################################
### Procedure to prompt for dark-pixel correction values
### for each band. Uses an auxillary dialog window to list
### input bands with a field for each to enter the value.
#######################################
proc SetDP () {
local string xmldp$, nirprtext$, redprtext$, labeltext$;
local numeric errdp;
### Create string variable with XML specification of dialog
### to enter dark-pixel-correction values
xmldp$ = '<?xml version="1.0"?>
<root>
<dialog id = "dpdlg" Title = "Set Dark Pixel Values" OnOK = "dpOK()" >
<label>Set Dark Pixel (Path Radiance)</label>
<label>Correction Values:</label>
<groupbox ExtraBorder = "3">
<pane id = "dplist" Orientation = "vertical"/>
</groupbox>
</dialog>
</root>';
### Parse XML text for dialog into memory; return error if there are syntax errors.
errdp = docdp.Parse(xmldp$);
if ( errdp < 0 ) {
PopupError( errdp ); # pop up an error dialog
Exit( );
}
#############################
### Modify the XML structure in
### memory before opening the dialog
##############################
### Get the id for the parent pane that will contain the list of bands and the
### numeric fields for the correction values
class XMLNODE dplist;
dplist = docdp.GetElementByID( "dplist" );
########################
### Add horizontal pane for the first row of dialog elements
### (label and numeric field for NIR band)
########################
class XMLNODE paneNIR;
paneNIR = dplist.NewChild( "pane" );
paneNIR.SetAttribute( "Orientation", "Horizontal");
# Add label with name of NIR band to the NIR pane
class XMLNODE labelNIR;
labelNIR = paneNIR.NewChild( "label" );
nirprtext$ = sprintf("%s.%s \ %s", FileNameGetName(nirfile$), FileNameGetExt(nirfile$), nirobjname$);
labelNIR.SetText(nirprtext$);
# Add editnumber field to the NIR pane and set its attributes
class XMLNODE prEditNIR;
prEditNIR = paneNIR.NewChild( "editnumber" );
prEditNIR.SetAttribute( "id", "prEditNIR" );
prEditNIR.SetAttribute( "Width", "5" );
prEditNIR.SetAttribute( "MinVal", "0" );
prEditNIR.SetAttribute( "MaxVal", "255" );
prEditNIR.SetAttribute( "Default", NumToStr(nirmin) );
prEditNIR.SetAttribute( "Precision", "0" );
########################
### Add horizontal pane for the second row of dialog elements
### (label and numeric field for RED band)
########################
class XMLNODE paneRED;
paneRED = dplist.NewChild( "pane" );
paneRED.SetAttribute( "Orientation", "Horizontal");
# Add label with name of RED band to the RED pane
class XMLNODE labelRED;
labelRED = paneRED.NewChild( "label" );
redprtext$ = sprintf("%s.%s \ %s", FileNameGetName(redfile$), FileNameGetExt(redfile$), redobjname$);
labelRED.SetText(redprtext$);
# Add editnumber field to the RED pane and set its attributes
class XMLNODE prEditRED;
prEditRED = paneRED.NewChild( "editnumber" );
prEditRED.SetAttribute( "id", "prEditRED" );
prEditRED.SetAttribute( "Width", "5" );
prEditRED.SetAttribute( "MinVal", "0" );
prEditRED.SetAttribute( "MaxVal", "255" );
prEditRED.SetAttribute( "Default", NumToStr(redmin) );
prEditRED.SetAttribute( "Precision", "0" );
#########################
### Loop to add rows for additional bands to be processed
#########################
if ( numbands > 0 ) { # prompt for values for additional bands
class XMLNODE panepr; # reusable XMLNODE class instances for new
class XMLNODE labelpr; # pane, label, and editnumber field for
class XMLNODE editpr; # each additional band
local string editid$; # string for id for editnumber dialog element
local numeric stringlistindex; # string number in stringlist (starts at 0)
for num = 1 to numbands {
stringlistindex = num - 1;
### create horizontal pane for row
panepr = dplist.NewChild( "pane" );
panepr.SetAttribute( "Orientation", "Horizontal" );
### get file/object names for current band (stored in stringlists)
infile$ = inrastfilelist.GetString( stringlistindex );
inobjname$ = inrastobjlist.GetString( stringlistindex );
labeltext$ = sprintf("%s.%s \ %s", FileNameGetName(infile$), FileNameGetExt(infile$), inobjname$);
### add label to band pane
labelpr = panepr.NewChild( "label" );
labelpr.SetText( labeltext$ );
### add editnumber field to the band pane and set its attributes
editpr = panepr.NewChild( "editnumber" );
editid$ = "editnum" + NumToStr( num );
editpr.SetAttribute( "id", editid$ );
editpr.SetAttribute( "Width", "5" );
editpr.SetAttribute( "MinVal", "0" );
editpr.SetAttribute( "MaxVal", "255" );
editpr.SetAttribute( "Default", NumToStr( rastmins[ num ] ) );
editpr.SetAttribute( "Precision", "0" );
}
}
### Get the dialog element from the parsed XML text and show error message if
### the dialog element can't be found
nodedp = docdp.GetElementByID( "dpdlg" );
if ( nodedp == 0 ) {
PopupMessage( "Could not find dialog node in XML document" );
Exit();
}
### Set the XMl dialog element as the source for the GUI_DLG class instance we
### are using for the dark-pixel-correction dialog window
dlgdp.SetXMLNode( nodedp );
### Open the dialog window as a modal dialog
dlgdp.DoModal();
} # end SetDP()
########################################
### Procedure to read dark-pixel correction values from
### dp dialog window when its OK button is pressed.
########################################
proc dpOK () {
dlgmain.GetCtrlByID("getfile").SetEnabled(1);
dlgmain.GetCtrlByID("getdir").SetEnabled(1);
local class GUI_FORMDATA dpdata;
dpdata = dlgdp.GetValues();
prNIR = dpdata.GetValueNum( "prEditNIR" );
# printf( "prNIR = %d\n", prNIR );
prRED = dpdata.GetValueNum( "prEditRED" );
# printf( "prRED = %d\n", prRED );
if ( numbands > 0 ) { # get values for additional bands and store in array
array numeric prband[ numbands ]; # array to hold dark pixel correction values for additional rasters
local string editid$; # string for id for editnumber dialog element
for num = 1 to numbands {
editid$ = "editnum" + NumToStr( num );
prband[ num ] = dpdata.GetValueNum( editid$ );
# printf("dp value for additional band %d = %d\n", num, prband[ num ] );
}
}
} # end dpOK ()
###################################
### Procedure to select output file
###################################
proc GetOutFile () {
infile$ = GetObjectFileName( NIR ); # directory path and input project file name
infilepath$ = FileNameGetPath( infile$ ); # gets path portion filename string
infilename$ = FileNameGetName( infile$ ); # gets filename portion of filename string
deffile$ = infilepath$ + "/" + infilename$ + "Out.rvc" ; # default output file name
outfile$ = GetOutputFileName( deffile$, "Select a Project File for script output:", "rvc" );
dlgmain.SetCtrlValueStr( "outfile", outfile$ );
} # end GetOutFile ()
#####################################
### Procedure to select directory in which TNTmips is installed
#####################################
proc GetMipsDir () {
mipspath$ = GetDirectory( "c:/Program Files/MicroImages/TNT67",
"Select directory in which TNTmips is installed:" );
cmapfilename$ = mipspath$ + "/" + "colormap.rvc";
cmapnum = SubObjectNumber( cmapfilename$, 0, "Rainbow1", "COLMAP" );
dlgmain.SetOkEnabled( 1 );
} # end GetMipsDir ()
#########################################
### Procedure to get the dialog window settings and store
### in a GUI_FORMDATA winsettings.
#########################################
proc GetValues () {
winsettings = dlgmain.GetValues();
devegNIR = winsettings.GetValueNum( "devegNIR" );
devegRED = winsettings.GetValueNum( "devegRED" );
numbands = winsettings.GetValueNum( "numbands" );
initavg = winsettings.GetValueNum( "avginitmed" );
smoothcurve = winsettings.GetValueNum( "smoothcurve" );
}
####################################################################
### Procedure to make a lookup table with the dark-pixel-corrected
### raster value for each possible input raster value. Corrected
### values less than 0 are set to 0. Lookup table is stored as a
### temporary raster with 1 line and 256 columns.
###############################################################
proc dkpixel ( raster BandRaw, raster Prlookup, numeric pr) {
# BandRaw = input imagery band to correct
# Prlookup = raster containing corrected band value for each raw value
# pr = value to subtract
local numeric i;
for i = 1 to 256 {
Prlookup[ 1, i ] = i - 1 - pr;
if ( Prlookup[ 1, i ] < 0 ) then
Prlookup[ 1, i ] = 0;
# printf("raster value %d = %d\n", i - 1, Prlookup[ 1, i ] );
}
} # end procedure dkpixel()
################################################################
### Procedure to compute and apply adjustments to band DN values
### for each vegetation index level.
################################################################
proc deveg ( raster Band, raster Prlookup, raster Corr2D, raster BandOut, CAD Graph, numeric initavg, numeric smoothcurve ) {
# Band = imagery band to devegetate
# Prlookup = temporary raster to use as lookup table for dark-pixel-corrected values
# Corr2D = saved correlation raster
# BandOut = output devegetated band raster
# Graph = CAD object graphing smoothed correlation curve used to make
# the band adjustment; for display over saved correlation raster
local numeric vegval, bandval, crastnum;
local raster TempCorr, BVavg;
#############################################################################
### tabulate for each vegetation index value the number of occurrences of
### each band DN value and store as 256 x 256 raster with line number = VI (0 - 255) + 1
### and column number = band DN value (0 to 255) + 1. Save correlation raster
### for later inspection.
#############################################################################
printf( " Creating correlation raster...\n" );
# Create temporary correlation raster with line number = VI level + 1
# and column number = band value + 1 (origin in upper left corner);
CreateTempRaster( TempCorr, 256, 256, "16-bit unsigned" );
for each VI {
vegval = VI;
bandval = Prlookup[ 1, Band + 1 ];
TempCorr[ vegval - 1, bandval - 1 ] += 1;
}
# make a saved version of the correlation raster rotated 90 degrees
# counterclockwise to put origin in lower left corner to match the CAD graph.
local numeric linout, colout;
local numeric linin, colin;
for linout = 1 to 256 {
for colout = 1 to 256 {
Corr2D[ linout, colout ] = TempCorr[ colout , 257 - linout ];
}
}
# copy standard color palette Rainbow1 to the output correlation raster Corr2D
crastnum = GetObjectNumber(Corr2D);
CopyObject( cmapfilename$, cmapnum, outfile$, crastnum );
class RVC_COLORMAP cmap;
class COLOR background;
background.red = 100; background.green = 100; background.blue = 100;
cmap = ColorMapFromRastVar(Corr2D);
ColorMapSetColor(cmap, 0, background);
ColorMapWriteToRastVar(Corr2D, cmap, "Rainbow1", "Rainbow1 color palette with 0 = white");
CreateHistogram(Corr2D);
CreatePyramid(Corr2D);
CloseRaster(Corr2D);
#########################################################################
#### Find smooth curve for plot of "average band value versus index level
#########################################################################
# Create temporary raster to hold the "average" band value for each index level
CreateTempRaster(BVavg, 1, 256, "16-bit unsigned");
printf( " Finding smooth curve through correlation values...\n\n" );
###############
# Compute initial "average" values for each index level as the median of
# average DN's within a 5-index spread centered on the current index level
###############
local numeric level; # VI level: (column number - 1) in BVavg raster and (line number - 1) in TempCorr
local numeric bv; # dn level in band = column number - 1 in TempCorr
local numeric bvavg; # average value of band value for 5-line window
local numeric total; # cumulative sum of average dn cell values
local numeric midpoint; # half of sum total
if ( initavg == 1 ) {
# Compute median for each VI by first finding for each bv the average count
# for five VI levels centered on the current level.
# Compute sum of all BV counts for each VI level; median is BV level where
# running sum of counts is half of total
for level = 1 to 256 {
total = 0;
midpoint = 0;
for bv = 1 to 256 {
bvavg = 0;
bvavg += TempCorr[ SetMax( 1, level-2 ), bv ];
bvavg += TempCorr[ SetMax( 1, level-1 ), bv ];
bvavg += TempCorr[ level, bv ];
bvavg += TempCorr[ SetMin( 256, level + 1 ), bv ];
bvavg += TempCorr[ SetMin( 226, level + 2 ), bv ];
bvavg /= 5;
total += bvavg;
}
midpoint = total / 2;
# printf("total = %d midpoint= %d\n", total, midpoint);
total = 0;
# loop again to find dn corresponding to midpoint
for bv = 1 to 256 {
bvavg = 0;
bvavg += TempCorr[ SetMax( 1, level-2 ), bv ];
bvavg += TempCorr[ SetMax( 1, level-1 ), bv ];
bvavg += TempCorr[ level, bv ];
bvavg += TempCorr[ SetMin( 256, level + 1 ), bv ];
bvavg += TempCorr[ SetMin( 226, level + 2 ), bv ];
bvavg /= 5;
total += bvavg;
if ( total >= midpoint ) {
BVavg[ 1, level ] = bv;
# printf( " Median for VI level %d = %d\n ", level, BVavg[ 1, level ] );
break;
}
}
}
}
else {
# Find median value for each single VI level with out averaging.
# Compute sum of all BV counts for each VI level; median is BV level where
# running sum of counts is half of total
for level = 0 to 255 {
total = 0;
midpoint = 0;
for bv = 1 to 256 {
total += TempCorr[ level + 1, bv ];
}
midpoint = total / 2;
#printf("total = %d midpoint= %d\n", total, midpoint);
total = 0;
# loop again to find dn corresponding to midpoint
for bv = 1 to 256 {
total += TempCorr[ level + 1, bv ];
if ( total >= midpoint ) {
BVavg[ 1, level + 1 ] = bv;
# printf( " Median for VI level %d = %d\n ", level, BVavg[ 1, level ] );
break;
}
}
}
}
if ( smoothcurve == 1) {
#####################
### Smooth curve by repeated focal-filtering of BVavg raster
#####################
### Create array to hold filter output
local array numeric filtout[256];
###########
### Sequence of progressively smaller median filters:
###########
# ### Median filter, filter width 23 (half-width 11)
# for level = 12 to 245 {
# filtout[ level ] = FocalMedian( BVavg, 1, 23, 1, level );
# }
#
# for level = 12 to 245 {
# BVavg[ 1, level ] = filtout[ level ]
# }
#
# ### Median filter, filter width 15 (half-width 7)
# for level = 8 to 249 {
# filtout[ level ] = FocalMedian( BVavg, 1, 15, 1, level );
# }
#
# for level = 8 to 249 {
# BVavg [ 1, level ] = filtout[ level ];
# }
#
# ### Median filter, filter width 9 (half-width 4)
# for level = 5 to 252 {
# filtout [ level ] = FocalMedian( BVavg, 1, 7, 1, level );
# }
#
# for level = 5 to 252 {
# BVavg [ 1, level ] = filtout[ level ];
# }
#
# ### Median filter, filter width 7 (half-width 3)
# for level = 4 to 253 {
# filtout [ level ] = FocalMedian( BVavg, 1, 7, 1, level );
# }
#
# for level = 4 to 253 {
# BVavg [ 1, level ] = filtout[ level ];
# }
### Median filter, filter width 5 (half-width 2)
for level = 3 to 254 {
filtout [ level ] = FocalMedian( BVavg, 1, 7, 1, level );
}
for level = 3 to 254 {
BVavg [ 1, level ] = filtout[ level ];
}
### Median filter, filter width 3 (half-width 1)
for level = 2 to 255 {
filtout [ level ] = FocalMedian( BVavg, 1, 3, 1, level );
}
for level = 2 to 255 {
BVavg [ 1, level ] = filtout [ level ];
}
####################
#### Sequence of progressively larger mean filters:
####################
### Mean filter, filter width = 3 (half-width = 1)
for level = 2 to 255 {
filtout[ level ] = FocalMean( BVavg, 1, 3, 1, level );
}
for level = 2 to 255 {
BVavg[ 1, level ] = filtout[ level ];
}
### Mean filter, filter width = 5 (half-width = 2)
for level = 3 to 254 {
filtout[ level ] = FocalMean( BVavg, 1, 5, 1, level );
}
for level = 3 to 254 {
BVavg[ 1, level ] = filtout[ level ];
}
}
for level = 0 to 255 {
printf( " Median for VI level %d = %d\n", level, BVavg[ 1, level + 1 ] );
}
###############################################################
####### Make CAD graph of band value versus VI level
###############################################################
# create arrays to hold x and y coordinates of graph line;
local array numeric graphx[256];
local array numeric graphy[256];
for level = 1 to 256 {
graphx[ level ] = level;
graphy[ level ] = BVavg[ 1, level ];
}
local class CADELEMOPT boxopt, graphopt;
local class COLOR boxcolor, graphcolor;
boxcolor.red = 20;
boxcolor.green = 20;
boxcolor.blue = 20;
graphcolor.red = 30;
graphcolor.green = 30;
graphcolor.blue = 30;
CADWriteBox( Graph, 1, 0, 0, 256, 256, 0, boxopt );
CADWriteLine( Graph, 1, 256, graphx, graphy, graphopt );
CloseCAD( Graph );
################################################################
### compute devegetation adjustment from DNlevel values
### and make output devegetated band raster
################################################################
local numeric bandmean; # global mean for band raster
local numeric band_deveg; # computed devegetated band value
bandmean = GlobalMean(Band);
printf( "\n computing devegetated band...\n\n" );
for each Band {
vegval = VI; # VI value for current Band raster cell position
bandval = Prlookup[ 1, Band + 1 ]; # dark-pixel-corrected band value from lookup raster
bvavg = BVavg[ 1, vegval + 1 ]; # "average" band value for current VI
BandOut = round( bandval * bandmean / bvavg ); # rescaled band value
}
SetNull(BandOut, 255);
CreateHistogram( BandOut );
CreatePyramid( BandOut );
CopySubobjects( Band, BandOut, "georef" );
CloseRaster( Band );
CloseRaster( BandOut );
CloseRaster( Corr2D );
CloseRaster( BVavg );
} ### end deveg()
######################################################################
######################################################################
#################### Main Program ##################################
######################################################################
######################################################################
clear();
#################################################
### Set up main dialog window
#################################################
### Create string variable with XML specification of main
### dialog enclosed in single quotes (for multiline text)
xmlmain$ = '<?xml version="1.0"?>
<root>
<dialog id = "maindlg" title="Suppress Vegetation by Forced Invariance" OnOpen = "OnOpenMain()" OnOK = "GetValues()" >
<label>Select input near-infrared and red wavelength bands:</label>
<pane Orientation = "Horizontal">
<pushbutton Name = "Select NIR..." OnPressed = "SelectNIR()" />
<edittext id = "NIRtext" Width = "40" ReadOnly = "True"/>
</pane>
<pane Orientation = "Horizontal">
<pushbutton id = "getred" Name = "Select RED..." Enabled = "0" OnPressed = "SelectRED()" />
<edittext id = "REDtext" Width = "40" ReadOnly = "True"/>
</pane>
<pane Orientation = "Horizontal">
<togglebutton id = "devegNIR" Name = "Devegetate NIR" HorizResize = "Fixed" />
<togglebutton id = "devegRED" Name = "Devegetate RED" HorizResize = "Fixed" />
</pane>
<pane Orientation = "Horizontal">
<editnumber id = "numbands" Default = "0" MinVal = "0" MaxVal = "10" Width = "3" HorizResize = "Fixed" Precision = "0" OnChanged = "OnNumChng()" />
<label>Number of additional bands to devegetate</label>
<pushbutton id = "getaddbands" Name = "Select Bands..." Enabled = "0" OnPressed = "GetAddBands()"/>
</pane>
<listbox id = "bandlist" Width = "40" Height = "5" Enabled = "0"/>
<pane Orientation = "Horizontal">
<pushbutton id = "setdp" Name = "Set dark pixel values..." VertResize = "Fixed" Enabled = "0" OnPressed = "SetDP()" />
<groupbox Name = " Curve-fitting Parameters: " ExtraBorder = "2">
<pane Orientation = "Vertical">
<togglebutton id = "avginitmed" Name = "Average initial medians"/>
<togglebutton id = "smoothcurve" Name = "Smooth Band vs VI curve"/>
</pane>
</groupbox>
</pane>
<pane Orientation = "Horizontal">
<pushbutton id = "getfile" Name = "Select Output File..." Enabled = "0" OnPressed = "GetOutFile()"/>
<edittext id = "outfile" Width = "50" ReadOnly = "True"/>
</pane>
<pane Orientation = "Horizontal">
<label HorizResize = "Fixed">Select directory in which TNTmips is installed:</label>
<pushbutton id = "getdir" Name = "Directory..." Enabled = "0" OnPressed = "GetMipsDir()"/>
</pane>
</dialog>
</root>';
### parse XML text for main dialog into memory; returns error code (number < 0)
### if there are syntax errors
errXML = docmain.Parse(xmlmain$);
if ( errXML < 0 ) {
PopupError( errXML ); # pop up an error dialog
Exit( );
}
### get the dialog element from the parsed XML text and show error
### message if the dialog element can't be found
nodemain = docmain.GetElementByID( "maindlg" );
if ( nodemain == 0 ) {
PopupMessage( "Could not find main dialog node in XML document" );
Exit();
}
### set the XML dialog element as the source for the GUI_DLG class instance
### we are using for the main dialog window
dlgmain.SetXMLNode(nodemain);
############################################
### Open main dialog window as modal dialog
############################################
dlgreturn = dlgmain.DoModal();
if ( dlgreturn == -1 ) {
Exit();
}
##################################
### Begin processing.
### Calculate vegetation index.
##################################
printf( "Computing vegetation index...\n\n" );
### Create temporary rasters for NIR and RED bands to use as lookup-tables for
### dark-pixel-corrected raster values
CreateTempRaster( NIRlookup, 1, 256, "8-bit unsigned" );
dkpixel(NIR, NIRlookup, prNIR );
CreateTempRaster( REDlookup, 1, 256, "8-bit unsigned" );
dkpixel(RED, REDlookup, prRED );
### create temporary raster for raw vegetation index
CreateTempRaster( RawVI, lins, cols, "32-bit float" );
### create output raster for rescaled integer vegetation index
CreateRaster( VI, outfile$, "VI", " ", lins, cols, "8-bit unsigned" );
SetNull( VI, 255 );
for each NIR {
# since we are using a lookup table (raster) to find the dark-pixel-corrected
# values to use in the NDVI calculation, we must check for
r = REDlookup[ 1, RED + 1 ];
n = NIRlookup[ 1, NIR + 1 ];
if ( IsNull( RED ) or IsNull( NIR ) or ( r == 0 and n == 0 ) ) then
RawVI = 340282346638528860000000000000000000000.000000;
else
RawVI = ( n - r ) / ( n + r );
}
SetNull(RawVI, 340282346638528860000000000000000000000.000000);
### Rescale RawVI to 8-bit unsigned integer range
### Use -0.25 as lower limit for rescaling. Set upper limit
### 10% of remaining range below the maximum
maxRawVI = GlobalMax( RawVI );
minRawVI = GlobalMin( RawVI );
printf( "Min raw NDVI = %0.6f\n", minRawVI );
printf( "Max raw NDVI = %0.6f\n", maxRawVI );
### compute scale factor and offset
lowlim = -0.25;
range = maxRawVI - lowlim;
uplim = maxRawVI - 0.1 * range;
printf( "computed upper limit raw NDVI = %0.6f\n", uplim );
printf( "preset lower limit raw NDVI = %0.6f\n", lowlim );
offset = abs( lowlim );
printf("offset = %0.6f\n\n", offset);
scale = 254 / ( uplim - lowlim );
printf( "scale = %0.5f\n\n", scale );
### rescale vegetation index to integer output raster VI
printf( "Rescaling vegetation index to unsigned integer range...\n\n" );
for each RawVI {
VI = floor( ( RawVI + offset ) * scale );
}
CloseRaster( RawVI );
if ( devegNIR == 0 ) {
CloseRaster( NIR );
CloseRaster( NIRlookup );
}
if ( devegRED == 0 ) {
CloseRaster( RED );
CloseRaster( REDlookup );
}
printf( " Max VI = %d\n", GlobalMax( VI ) );
printf( " Min VI = %d\n\n", GlobalMin( VI ) );
CreatePyramid( VI );
CreateHistogram( VI );
CopySubobjects( NIR, VI, "georef" );
############ Compute devegetation adjustments for each band
############ and make new devegetated band raster in output file.
if ( devegNIR == 1 ) {
printf( "Devegetating NIR Band...\n" );
nirdvobjname$ = nirobjname$ + "DV";
CreateRaster( Corr2Dnir, outfile$, "Corr2Dnir", "2D correlation VI vs NIR", 256, 256, "16-bit unsigned" );
CreateRaster( NIRdeveg, outfile$, nirdvobjname$, "Devegetated NIR", lins, cols, datatype$ );
CreateCAD( GraphNIR, outfile$, "GraphVIvNIR", "CAD graph of smoothed Band vs Index" );
deveg( NIR, NIRlookup, Corr2Dnir, NIRdeveg, GraphNIR, initavg, smoothcurve );
}
if ( devegRED == 1 ) {
printf( "Devegetating Red Band...\n" );
reddvobjname$ = redobjname$ + "DV";
CreateRaster( Corr2Dred, outfile$, "Corr2Dred", "2D correlation VI vs Red", 256, 256, "16-bit unsigned" );
CreateRaster( REDdeveg, outfile$, reddvobjname$, "Devegetated Red", lins, cols, datatype$ );
CreateCAD( GraphRED, outfile$, "GraphVIvRED", "CAD graph of smoothed Band vs Index" );
deveg( RED, REDlookup, Corr2Dred, REDdeveg, GraphRED, initavg, smoothcurve );
}
if ( numbands > 0 ) {
### Create a temporary raster for band to use as lookup-table for
### dark-pixel corrected raster values
CreateTempRaster( PRlookup, 1, 256, "8-bit unsigned" );
for num = 1 to numbands {
num$ = NumToStr( num );
printf( "Devegetating additional band %s...\n ", num$ );
infile$ = inrastfilelist.GetString( num - 1 );
inobjname$ = inrastobjlist.GetString( num - 1 );
OpenRaster( Rin, infile$, inobjname$ );
### create correlation raster
corr$ = "Corr2D" + inobjname$;
corrdesc$ = "2D correlation VI vs " + inobjname$;
CreateRaster( Corr2D, outfile$, corr$, corrdesc$, 256, 256, "16-bit unsigned" );
### create raster for devegetated band
outobjname$ = inobjname$ + "DV";
CreateRaster( Deveg, outfile$, outobjname$, "Devegetated band", lins, cols, datatype$ );
cadname$ = "GraphVIv" + inobjname$;
CreateCAD( Graph, outfile$, cadname$, "CAD graph of Band vs Index" );
### call procedure to make dark-pixel-correction lookup table
dkpixel(Rin, PRlookup, prband[ num ] );
### call procedure to devegetate the band
deveg(Rin, PRlookup, Corr2D, Deveg, Graph, initavg, smoothcurve );
}
}
CloseRaster( VI );
printf( "Done" );
Download