#### 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 class RVC_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" );