# PipelineNDVIfromTIFF.sml

# Sample standalone script to demonstrate an image pipeline application.

# Computes a Normalized Difference Vegetation Index (NDVI) raster from Red and
# Near-infrared bands directly from an input GeoTIFF file containing a 4-band multispectral
# QuickBird or Ikonos satellite image.  The GeoTIFF file does not need to have been previously linked
# or imported to TNTmips.  

# NDVI is calculated for each image cell using the formula:
# 				1000 * (Near Infrared - Red) / (Near Infrared + Red)
# Values in the resulting NDVI raster can range from -1000 to + 1000 (using Signed 16-bit raster type).

# If the input image is projected, there may be 0-value cells around the edges.  A ValidityNear filter
# is used to identify cells that have zero value in both bands; the filter creates a null mask
# for each band with such cells marked as null.

# Path radiance (atmospheric scattering / haze) is an additive effect that varies between
# bands, and so should be removed before calculating NDVI.  This script uses the minimum (non-null)
# raster value in each band as an approximation to the Path Radiance value for that band. This value
# is subtracted using the SCALEOFFSET filter (scale = 1.00, offset value negative).

# The LINEAR filter (linear combination of band values) is used to perform the addition and
# subtraction of the band values.

# The NDVI image is produced as a single 16-bit signed integer RVC raster in a Project File.     

# 8 January 2008
# Randy Smith, MicroImages, Inc.

# Requires Version 2008:74 or later of the TNT products.

numeric err;

# error checking procedure
proc ReportError(numeric linenum, numeric err) {
	printf("FAILED -line: %d, error: %d\n", linenum - 1, err);
	PopupError(err); 
	}


################################################################################
#################################  Main program  ###############################
clear();

# set context so script does not exit on error, allowing manual
# handling of errors.
_context.AbortOnError = 0;


# CHOOSE GEOTIFF FILE containing QuickBird or IKONOS 4-band multispectral image
class FILEPATH inFilePath(GetInputFileName("","Choose GeoTIFF file with source image:", ".tif"));

# CHOOSE OUTPUT RASTER
# DlgGetObject creates only an ObjItem for the output; the new raster will then be created
# by the TARGET_RVC with the desired compression.  The ObjectPath of the ObjItem is not
# valid until the output raster is actually created by the pipeline.
class RVC_OBJITEM ndviOutObjItem;
DlgGetObject("Choose new raster for NDVI result:", "Raster", ndviOutObjItem, "NewOnly");

# PIPELINE SOURCE
# set GeoTIFF file containing a four-band QuickBird or Ikonos image as source
class IMAGE_PIPELINE_SOURCE_TIFF source_tif( inFilePath );
err = source_tif.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Source TIFF file initialized");
	}

if (source_tif.GetNumSamples() == 4) 
	printf("Number of bands in TIFF file = %d\n", source_tif.GetNumSamples() );
else {
	print("Selected TIFF file does not contain 4 bands; exiting now.");
	Exit(); 
	}

# Get georeference from the GeoTIFF file to print what CRS is used.
class IMAGE_PIPELINE_GEOREFERENCE georef;
georef = source_tif.GetGeoreference();

class SR_COORDREFSYS crs = georef.GetCRS();
printf("Coordinate reference system = %s\n", crs.Name );


# PIPELINE FILTER_SELECT to select red and near-infrared bands for processing (don't need to carry the rest
# forward through the pipeline).
# The source TIFF file consists of 4 bands (samples) per pixel, which have numeric indices 0 = blue,
# 1 = green, 2 = red, 3 = near-infrared.  These samples are selected by index number in this filter.
numeric numSamples = 2;		# number of bands from 4-band TIFF file to be used
array numeric srcSamples[2];	# array to contain the numeric indices of the bands to select
srcSamples[1] = 2;		# index for red band
srcSamples[2] = 3;		# index for near-infrared band

class IMAGE_PIPELINE_FILTER_SELECT filter_select(source_tif, srcSamples, numSamples);
err = filter_select.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized band selection filter");
	}

# the result of this FILTER_SELECT is an "image" with two samples per pixel, with the red value as index 0
# and the near-infrared value as index 1.

class STRING pixelType = filter_select.GetPixelType();
printf("Pixel type in image = %s\n", pixelType);

class STRING datatypeRed$, datatypeNIR$;
class IMAGE_PIPELINE_PIXEL pixel;
pixel = filter_select.GetPixelValueMax();
datatypeRed$ = pixel.GetSample(0).GetDataType();
datatypeNIR$ = pixel.GetSample(1).GetDataType();
printf("Datatype of red band = %s\n", datatypeRed$ );
printf("Datatype of near-infrared band = %s\n", datatypeNIR$ );



# PIPELINE FILTER_VALIDITYNEAR
# Filter to set null for cells based on nearness (distance) to a specified value
# In this case test for zero value in both bands (assumed to be non-image cells) and set
# such cells to null (0 value in null mask created by the filter)

# validity flag 
numeric resultValid = 0;  # set 1 if pixels near test value should be valid, 0 if invalid;
									# we want 0 values to be invalid, so set this parameter to false (0)
numeric distance = 0.5;		# distance to test must be greater than 0; use a decimal value less
									# than 1 to test for exact match to an integer value
class IMAGE_PIPELINE_FILTER_VALIDITYNEAR validity(filter_select, "Orthogonal", resultValid); 

# define two 0-value samples to provide the set of values 
class IMAGE_PIPELINE_SAMPLE redNullSamp("Gray", "UINT16", 0);
class IMAGE_PIPELINE_SAMPLE nirNullSamp("Gray", "UINT16", 0);
# define a text pixel with the two 0-value samples to pass to the validity filter
class IMAGE_PIPELINE_PIXEL testPixel(redNullSamp, nirNullSamp);
# add the pixel value to test and the test distance to the validity filter;
validity.AddValue(testPixel, distance);  

array numeric samplesToTest[2];	# array to set which image samples to test (in this case, both)
samplesToTest[1] = 0;	# sample index for red
samplesToTest[2] = 1;	# sample index for near-infrared
validity.SetSamplesTest(samplesToTest, 2);

err = validity.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized ValidityNear filter.");
	}


# GET MINIMUM AND MAXIMUM VALUES of the red and near-infrared bands.
# First call method on stage to compute the range of values in the bands; using the "Exact"
# value for the RangeType parameter in this method to compute the actual range of values using the 
# rasters without sampling (i.e. computing using the pyramids, if any).  The "Default" value returns
# the range of the raster datatype.
validity.ComputePixelRange("Exact"); 
class IMAGE_PIPELINE_PIXEL minVals, maxVals;		# container to hold the minimum and maximum value for each band
minVals = validity.GetPixelValueMin();
maxVals = validity.GetPixelValueMax(); 

printf("Minimum value of red band = %d\n", minVals.GetSample(0).GetValue() );
printf("Maximum value of red band = %d\n", maxVals.GetSample(0).GetValue() );

printf("Minimum value of near-infrared band = %d\n", minVals.GetSample(1).GetValue() );
printf("Maximum value of near-infrared band = %d\n", maxVals.GetSample(1).GetValue() );



# PIPELINE FILTER_SCALEOFFSET applies scale (multiply by constant) and offset (add constant value) to each cell
# Use with scale = 1 and negative offset to subtract the minimum valid value from each of the bands
# as an approximate correction for path-radiance.
class IMAGE_PIPELINE_FILTER_SCALEOFFSET offsetMin(validity);

# set scale = 1 and negative offset equal to minimum value for red band (0) and near-infrared band (1)
offsetMin.SetScaleOffset(0, 1, minVals.GetSample(0).GetValue() * -1);
offsetMin.SetScaleOffset(1, 1, minVals.GetSample(1).GetValue() * -1);
err = offsetMin.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized Scale-offset filter for path-radiance correction.");
	}
 

# PIPELINE FILTER_LINEAR produces a linear combination of the band values using a coefficient for each band
# Cell Value = a * Band1value  +  b * Band2value + ...
# Use here to subtract band values to provide dividend for the NDVI expression: NIR - Red; 
# Set number of samples in target to 1; coefficient for NIR = 1 and for Red = -1
numeric numTargetSamples = 1;
numeric targetSampleNum = 0;	# 0th target sample
numeric sourceSampleNum = 0;	# 0th source sample (red)
numeric coefValue;		# value of the coefficient

class IMAGE_PIPELINE_FILTER_LINEAR linearSubtract(offsetMin, numTargetSamples);
coefValue = -1;	# set coefficient value = -1 for red 
linearSubtract.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue); 
sourceSampleNum = 1;	# 1st (near-infrared) source sample
coefValue = 1;		# set coefficient value = 1 for near-infrared
linearSubtract.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue);
err = linearSubtract.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized linear filter for subtracting red from near-infrared.");
	}	 


# PIPELINE FILTER_LINEAR
# Use here to add band values to provide divisor for the NDVI expression: NIR + Red;
# Set number of samples in target to 1; coefficient for both NIR and Red = 1;
class IMAGE_PIPELINE_FILTER_LINEAR linearAdd(offsetMin, numTargetSamples);
sourceSampleNum = 0;	# 0th source sample (red)
linearAdd.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue);
sourceSampleNum = 1;	# 1st source sample (near-infrared)
linearAdd.SetCoefficient(targetSampleNum, sourceSampleNum, coefValue);
err = linearAdd.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized linear filter for adding red and near-infrared.");
	}


# PIPELINE FILTER_DIVIDE: filter to divide two images
# Use to perform the final NDVI division using the results from the two Linear filters
# The stage takes two source stages, the dividend stage followed by the divisor stage,
# Output of this filter is a single 64-bit floating-point image with possible range -1.00 to + 1.00
class IMAGE_PIPELINE_FILTER_DIVIDE divide(linearSubtract, linearAdd);
err = divide.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized divide filter to compute NDVI.");
	}


# PIPELINE FILTER_SCALEOFFSET
# Scale floating-point values to range -1000 to + 1000 using scale factor = 1000 and offset = 0
class IMAGE_PIPELINE_FILTER_SCALEOFFSET scale1000(divide);
scale1000.SetScaleOffset(0, 1000, 0);
err = scale1000.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized scale-offset filter to rescale NDVI.");
	}



# PIPELINE FILTER_DATATYPE: filter to change datatype
# Change rescaled NDVI from 64-bit floating point to 16-bit signed integer
class IMAGE_PIPELINE_FILTER_DATATYPE changeToSINT16(scale1000, "SINT16", "Invalidate");
err = changeToSINT16.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized datatype filter to convert to signed 16-bit.");
	}


# PIPELINE TARGET
class IMAGE_PIPELINE_TARGET_RVC target_rvc(changeToSINT16, ndviOutObjItem);
target_rvc.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized NDVI target.");
	}


# EXECUTE the pipeline process
print("Processing...");
target_rvc.Process();

print("Done.");