# PipelineNDVIfromTIFF_CropAndMaskFromRegion.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 computations are restricted to areas of the image delineated by
# polygons in an input vector object.  The polygons are converted to a region which is used to
# mask the areas outside the polygons and to crop the image to the extent of the polygons.  

# 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.     

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



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

printf("Number of bands in TIFF file = %d\n", source_tif.GetNumSamples() );

if (source_tif.GetNumSamples() <> 4) {
	print("Selected TIFF file does not contain 4 bands; exiting now.");
	Exit(); 
	}

# Get georeference and coordinate reference system from the GeoTIFF image.
class IMAGE_PIPELINE_GEOREFERENCE georefImg;
georefImg = source_tif.GetGeoreference();

class SR_COORDREFSYS imgCRS = georefImg.GetCRS();

if (!imgCRS.IsDefined() ) {
	print("Image is not georeferenced; exiting now.");
	Exit();
	}
else if (imgCRS.IsLocal() ) {
	print("Image has local engineering georeference; exiting now.");
	Exit();
	}
else
	printf("Image coordinate reference system = %s\n", imgCRS.Name );

# Get coordinate transformation from image to its defined CRS
class TRANS2D_MAPGEN imgToCRS = georefImg.GetTransGen();



# CHOOSE VECTOR OBJECT WITH POLYGONS TO INDICATE AREAS TO PROCESS
class RVC_VECTOR polyVect;
class RVC_OBJITEM vectObjItem;
class RVC_GEOREFERENCE vectGeoref;
class SR_COORDREFSYS vectCRS;
DlgGetObject("Choose vector object with polygons outlining areas to process:", "Vector", vectObjItem, "ExistingOnly");
polyVect.Open(vectObjItem, "Read");

if (polyVect.$Info.NumPolys == 0) {
	print("Vector object contains no polygons; exiting now.");
	Exit();
	}

err = polyVect.GetDefaultGeoref(vectGeoref);
if (err < 0 ) {
	ReportError(_context.CurrentLineNum, err);
	print("Polygon vector must be georeferenced and is not; exiting now.");
	Exit();
	}

vectCRS = vectGeoref.GetCoordRefSys();
printf("Vector coordinate reference system = %s\n", vectCRS.Name);



# 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 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$ );



# CONVERT INPUT VECTOR TO A REGION FOR USE AS A PIPELINE SOURCE
class REGION2D RegFromVect = ConvertVectToRegion(polyVect, GetLastUsedGeorefObject(polyVect) );

# if vector has different coordinate reference system from image, reproject region
if (vectCRS.Name <> imgCRS.Name) {
	printf("Converting region to image CRS: %s\n", imgCRS.Name);
	RegFromVect.ConvertTo(imgCRS);
	}

class RECT vecRegionExtents = RegFromVect.Extents;
printf("Minimum map extents of polyVect region: %.2f, %.2f\n", vecRegionExtents.x1, vecRegionExtents.y1);
printf("Maximum map extents of polyVect region: %.2f, %.2f\n\n", vecRegionExtents.x2, vecRegionExtents.y2);



# GET EXTENTS OF THE IMAGE SOURCE AS A REGION
class REGION2D imgReg;
err = filter_select.ComputeGeoreferenceRegion(imgReg);
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);



# CHECK THAT REGION FROM VECTOR IS CONTAINED WITHIN IMAGE REGION
if (imgReg.TestRegion(RegFromVect, "FullInside") ) {
	print("Image region contains vector extents.");
	}
else {
	print("Input vector is not contained within the source image.  Exiting now.");
	Exit();
	}



# PIPELINE SOURCE FOR REGION CREATED FROM THE VECTOR POLYGONS
# Construct using TRANS2D_MAPGEN with coordinate transformation from source image to its CRS.
# Use source image as a reference so the "image dimensions" of the region source match,
# enabling the region source to be used as a mask in the MASKVALIDITY filter
class IMAGE_PIPELINE_SOURCE_REGION sourceReg(RegFromVect, imgToCRS, filter_select);
err = sourceReg.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized region source.");
	printf("Size of region = %d lines, %d columns\n", sourceReg.GetTotalRows(), sourceReg.GetTotalColumns() );
	}

# Get the extents of the region as a RECT in image coordinates
# to use to crop the source image.  These are the extents of the region used to
# create the region source, not those of the image used as its reference.
class RECT regExtents;
regExtents = sourceReg.GetRegionExtents();
printf("Minimum values of cropping rectangle from polygons: x = %.2f, y = %.2f\n", regExtents.x1, regExtents.y1);
printf("Maximum values of cropping rectangle from polygons: x = %.2f, y = %.2f\n", regExtents.x2, regExtents.y2);



# SET UP PIPELINE FILTER TO MASK PORTIONS OF SOURCE IMAGE
# mask area outside the region created from the vector polygons
class IMAGE_PIPELINE_FILTER_MASKVALIDITY filterMask(filter_select, sourceReg);
err = filterMask.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized image mask filter.");
	}


# SET UP PIPELINE FILTER TO CROP THE SOURCE IMAGE
class IMAGE_PIPELINE_FILTER_CROP	filter_Crop(filterMask, regExtents);
err = filter_Crop.Initialize();
if (err < 0 )
	ReportError(_context.CurrentLineNum, err);
else {
	print("Initialized image crop filter.");
	}



# 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_Crop, "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.");