# REGSTATp.SML - Allows user to select an area of a non-composite raster
# by selecting and deselecting polygons from a vector layer over the top of it.
# The script then outputs the number of cells, number of non-null cells,
# minimum, maximum, mean, standard deviation, area, perimeter, centroid,
# and surface area over that region. The user may choose any distance or
# area units for the output to be displayed in. 
#
# Requires TNTmips version 7.0
# Updated 29 April 2005


# The following symbols are predefined
#    class VIEW View            {use to access the view the tool script is attached to}
#    class GROUP Group          {use to access the group being viewed if the script is run from a group view}

# Variable declarations
class XmForm form, buttonRow, plusMinus;
class XmSeparator line1, line2;
class GRE_LAYER_VECTOR vectorLayer;
class GRE_LAYER rasterLayer;
class Vector targetVector;
class Raster targetRaster;
class XmDrawingArea da;
class GC gc;
class PointTool pointTool;
string vectorName$, rasterName$, regionMode$;
class PushButtonItem saveButton, closeButton;
class XmPushButton clearicon;
class ToggleButtonItem plusicon,  minusicon;
class XmOptionMenu distMenu, areaMenu;
class REGION2D reg;
numeric min, max, mean, stdDev, count, cells, sum, sumsqr, surface, area, perimeter;
class POINT2D centroid;
numeric setDefaultWhenClose = false;
numeric distScale, areaScale, numPolys;
array numeric polyarray[1000];
class GRE_GROUP group;

# Checks layer to see if it is valid.
func checkLayer() {
	local numeric valid = false;

	# Get names layers if usable.  If not output error messages.
	if (group.ActiveLayer.TypeID == "") {
		rasterName$ = "Group has no layers!";
		}
	else {
		if (group.LastLayer.TypeID == "Vector") {
			vectorLayer = group.LastLayer;
			DispGetVectorFromLayer(targetVector, vectorLayer);
			if (targetVector.$Info.NumPolys > 0) {
				vectorName$ = vectorLayer.Name;
				valid = true;
				}
			else vectorName$ = "No polygons!";
			}
		else vectorName$ = "Not a vector!";
		if (group.FirstLayer.TypeID == "Raster") {
			rasterLayer = group.FirstLayer;
			DispGetRasterFromLayer(targetRaster, rasterLayer);
			if (targetRaster.$Info.Type != "binary" and
				targetRaster.$Info.Type != "4-bit unsigned" and
				targetRaster.$Info.Type != "8-bit signed" and
				targetRaster.$Info.Type != "8-bit unsigned" and
				targetRaster.$Info.Type != "16-bit signed" and
				targetRaster.$Info.Type != "16-bit unsigned" and
				targetRaster.$Info.Type != "32-bit signed" and
				targetRaster.$Info.Type != "32-bit unsigned" and
				targetRaster.$Info.Type != "32-bit floating-point" and
				targetRaster.$Info.Type != "64-bit signed" and
				targetRaster.$Info.Type != "64-bit unsigned" and
				targetRaster.$Info.Type != "64-bit floating-point") {
				rasterName$ = "Type not supported!";
				valid = false;
				}
			else
				rasterName$ = rasterLayer.Name;
			}
		else {
			rasterName$ = "Not a raster!";
			valid = false;
			}
		}
	return valid;
	}

# Callback for drawing area expose.  Draws text.
proc cbRedraw() {
	local numeric larea, lperimeter, lsurface;

	# Scale the values to specified units.
	larea = areaScale * area;
	lperimeter = distScale * perimeter;
	lsurface = areaScale * surface;
	if (gc == 0) return;
	ActivateGC(gc);

	# Clear the drawing area and redraw text.
	SetColorName("gray75");
	FillRect(0, 0, da.width, da.height);
	SetColorName("black");
	if (cells > 0) {
		DrawInterfaceText(sprintf("Vector: %s\nRaster: %s\nCells: %d\nNull Cells: %d\nMinimum: %.2f\nMaximum: %.2f\nMean: %.2f\nStandard Deviation: %.2f\nArea: %.2f\nPerimeter: %.2f\nCentroid: %.2f, %.2f\nSurface Area: %.2f",vectorName$, rasterName$, count, cells - count, min, max, mean, stdDev, larea, lperimeter, centroid.x, centroid.y, lsurface), 0, 10);
		}
	else DrawInterfaceText(sprintf("Vector: %s\nRaster: %s\nCells:\nNull Cells:\nMinimum:\nMaximum:\nMean:\nStandard Deviation:\nArea:\nPerimeter:\nCentroid:\nSurface Area:",vectorName$, rasterName$), 0, 10);
	}

# Callback for when the active layer changes.
proc cbLayer() {
	checkLayer();
	cbRedraw();
	}

# Callback for when the active group changes.
proc cbGroup() {
	group = Layout.ActiveGroup;
	WidgetAddCallback(group.LayerSelectedCallback, cbLayer);
	cbLayer();
	}

# Callback for when user clicks the right mouse button on the polygon tool
# (or clicks apply).
proc cbToolApply(class pointTool tool) {
	# If the selected layer is not valid, don't do anything.
	if (checkLayer()) {
		# Set local variables
		local numeric sum, sumsqr, xscale, yscale, zscale;
		local numeric current, right, down, downright, elemNum;
		local class POINT2D point;
		local class StatusHandle status;
		local class StatusContext context;
		cells = 0; min = 0; max = 0; mean = 0; stdDev = 0; sum = 0; sumsqr = 0; count = 0; surface = 0; area = 0; perimeter = 0;
		centroid.x = 0; centroid.y = 0; current = 0; right = 0; down = 0; downright = 0;
		xscale = ColScale(targetRaster);
		yscale = LinScale(targetRaster);
		zscale = group.FirstLayer.zscale;

		# Check point.
		point.x = tool.Point.x;
		point.y = tool.Point.y;

		point = TransPoint2D(point, ViewGetTransViewToScreen(View, 1));
		point = TransPoint2D(point, ViewGetTransMapToView(View, vectorLayer.Projection, 1));

		numeric elementNum = FindClosestPoly(targetVector, point.x, point.y, GetLastUsedGeorefObject(targetVector));

		if (elementNum > 0) {
			if (regionMode$ == "plus") {
				# Add to polygon list.
				numPolys += 1;
				polyarray[numPolys] = elementNum;
				vectorLayer.Poly.HighlightSingle(elementNum, "Add");
				}
			else {
				if (numPolys > 1) {
					# Remove from polygon list.
					local numeric i;
					for i = 1 to numPolys
						if (polyarray[i] == elementNum) {
							polyarray[i] = polyarray[numPolys];
							numPolys -= 1;
							vectorLayer.Poly.HighlightSingle(elementNum, "Subtract");
							break;
							}
					}
				else if (numPolys == 1 and polyarray[1] == elementNum) {
					vectorLayer.Poly.HighlightSingle(elementNum, "Subtract");
					numPolys = 0;
					}
				}

			# Get region and transform it to the appropriate coordinate system.
			if (numPolys > 0) {
				reg = ConvertVectorPolysToRegion(targetVector, GetLastUsedGeorefObject(targetVector), polyarray, numPolys);
				}
			else ClearRegion(reg);

			# Create a status bar so user knows the script is doing something.
			status = StatusDialogCreate(form);
			context = StatusContextCreate(status);
			StatusSetMessage(context, "Computing values...");

			# Loop on region to calculate statistics.
			local numeric lin, col;
			foreach targetRaster[lin, col] in reg {
				if (!IsNull(targetRaster)) {
					if (count == 0) {
						max = targetRaster;
						min = targetRaster;
						}
					else if (targetRaster > max) {
						max = targetRaster;
						}
					else if (targetRaster < min) {
						min = targetRaster;
						}
					sum += targetRaster;
					sumsqr += sqr(targetRaster);
					count += 1;
					}

				# Surface area is estimated by by adding the areas of triangles created
				# by the current cell, the cell below it, the cell to the right of it,
				# and the cell to the lower right of it.
				# Estimate over null cells by keeping the last known value.
				if (!IsNull(targetRaster))
					current = targetRaster;
				if (!IsNull(targetRaster[lin,col+1]))
					right = targetRaster[lin,col+1];
				if (!IsNull(targetRaster[lin+1,col]))
					down = targetRaster[lin+1,col];
				if (!IsNull(targetRaster[lin+1,col+1]))
					downright = targetRaster[lin+1,col+1];
				surface += .5*sqrt(sqr(yscale*current*zscale-yscale*right*zscale)
								+sqr(xscale*current*zscale-xscale*down*zscale)
								+sqr(xscale*yscale));
				surface += .5*sqrt(sqr(yscale*downright*zscale-yscale*right*zscale)
								+sqr(xscale*downright*zscale-xscale*down*zscale)
								+sqr(xscale*yscale));
				cells += 1;
				}
			CloseRaster(targetRaster);

			# Avoid division by zero when calculating mean and standard deviation.
			if (count > 1) {
				mean = sum / count;
				stdDev = sqrt((sumsqr - sqr(sum) / count) / (count - 1));
				area = reg.GetArea();
				perimeter = reg.GetPerimeter();
				centroid = reg.GetCentroid();
				}

			cbRedraw();

			# Destroy the status bar
			StatusContextDestroy(context);
			StatusDialogDestroy(status);
			}
		}
	}

# Called when a unit type is selected.
proc cbDistUnits(class Widget parent, unit$) {
	distScale = GetUnitConvDist("meters", unit$);
	cbRedraw();
	}

# Called when a unit type is selected.
proc cbAreaUnits(class Widget parent, unit$) {
	areaScale = GetUnitConvArea("square meters", unit$);
	cbRedraw();
	}

# Called when the save button is pressed.  Saves statistics to text file.
proc cbSave() {
	local class FILE outfile;
	local numeric larea, lperimeter, lsurface;

	# Scale the values to specified units.
	larea = areaScale * area;
	lperimeter = distScale * perimeter;
	lsurface = areaScale * surface;

	# Open output file, write to it, and close it
	outfile = GetOutputTextFile(_context.ScriptDir, "Choose file...", "txt");
	fprintf(outfile, "Raster: %s\nCells: %d\nNull Cells: %d\nMinimum: %.2f\nMaximum: %.2f\nMean: %.2f\nStandard Deviation: %.2f\nArea: %.2f\nPerimeter: %.2f\nCentroid: %.2f, %.2f\nSurface Area: %.2f",
							rasterName$, count, cells - count, min, max, mean, stdDev, larea, lperimeter, centroid.x, centroid.y, lsurface);
	fprintf(outfile, "\nDistance Units: %s\nArea Units: %s\n\n", distMenu.value, areaMenu.value);
	fclose(outfile);
	}

# Called when the user presses the clear region button.
proc cbClearRegion() {
	vectorLayer.UnHighlightAllElements();
	numPolys = 0;
	}

# Called when the close button is pressed.  Closes the dialogs.
proc cbClose() {
	pointTool.Managed = 0;
	pointTool.HasPosition = 0;
	DialogClose(form);
	if (setDefaultWhenClose) {
		setDefaultWhenClose = false;
		View.SetDefaultTool();
		}
	}

# Called when the user selects plus region.
proc cbPlusRegion() {
	regionMode$ = "plus";
	}

# Called when the user selects minus region.
proc cbMinusRegion() {
	regionMode$ = "minus";
	}

# Called the first time the tool is activated.
# If the tool implements a dialog it should be created (but not displayed) here.
func OnInitialize () {
	if (Layout) {
		WidgetAddCallback(Layout.GroupSelectedCallback, cbGroup);
		group = Layout.ActiveGroup;
		}
	else
		group = Group;
	WidgetAddCallback(group.LayerSelectedCallback, cbLayer);

	# Set up tool dialog.
	form = CreateFormDialog("Region Statistics");
	form.marginHeight = 2;
	form.marginWidth = 2;
	WidgetAddCallback(form.Shell.PopdownCallback, cbClose);

	# add region mode (+,-) region icons
	# set the icons to initial state - plus region
	plusicon = CreateToggleButtonItem("Add to region", cbPlusRegion);
	plusicon.IconClass = "standard";
	plusicon.IconName = "add_select";
	plusicon.OneOfMany = true;

	minusicon = CreateToggleButtonItem("Subtract from region", cbMinusRegion);
	minusicon.IconClass = "standard";
	minusicon.IconName = "minus_select";
	minusicon.OneOfMany = true;

	plusMinus = CreateForm(form);
	CreateIconButtonRow(plusMinus, plusicon, minusicon);
	plusMinus.TopWidget = form;
	plusMinus.LeftWidget = form;

	clearicon = CreateIconPushButton(form, "standard", "delete", "Clear regions");
	clearicon.TopWidget = form;
	clearicon.leftWidget = plusMinus;
	WidgetAddCallback(clearicon.ActivateCallback, cbClearRegion);

	# Drawing area for displaying text.
	da = CreateDrawingArea(form, 185, 301);
	da.topWidget = plusMinus;
	da.leftWidget = form;
	da.rightWidget = form;
	WidgetAddCallback(da.ExposeCallback, cbRedraw);

	# Separator between output and unit tools.
	line1 = CreateHorizontalSeparator(form);
	line1.topWidget = da;
	line1.leftWidget = form;
	line1.rightWidget = form;
	line1.bottomOffset = 2;

	# Menus for selecting units.
	distMenu = CreateUnitOptionMenu(form, "distance_units_c", cbDistUnits, 2, 0);
	distMenu.topWidget = line1;
	distMenu.leftWidget = form;

	areaMenu = CreateUnitOptionMenu(form, "area_units_c", cbAreaUnits, 1, 0);
	areaMenu.topWidget = distMenu;
	areaMenu.leftWidget = form;

	# Separator between unit tools and buttons.
	line2 = CreateHorizontalSeparator(form);
	line2.topWidget = areaMenu;
	line2.leftWidget = form;
	line2.rightWidget = form;
	line2.topOffset = 2;

	# Allows user to save statistics to text file.
	saveButton = CreatePushButtonItem("Save As...", cbSave);

	# Closes the dialog.
	closeButton = CreatePushButtonItem("Close", cbClose);

	# Row of buttons.
	buttonRow = CreateButtonRow(form, saveButton, closeButton);
	buttonRow.topWidget = line2;
	buttonRow.leftWidget = form;
	buttonRow.rightWidget = form;

	# Add point tool.
	pointTool = ViewCreatePointTool(View);
	ToolAddCallback(pointTool.ActivateCallback, cbToolApply);

	}  # end of OnInitialize 
 
# Called when tool is to be destroyed, will not be called if tool was never activated.
func OnDestroy () {
	pointTool.Managed = 0;
	DestroyGC(gc);
	DestroyWidget(form);
	}  # end of OnDestroy 
 
# Called when tool is activated.
func OnActivate () {
	pointTool.Managed = 1;
	pointTool.HasPosition = 0;
	plusicon.Selected = 1;
	if (checkLayer())
		cbClearRegion();
	DialogOpen(form);
	if (gc == 0)
		gc = CreateGCForDrawingArea(da);
	cells = 0; min = 0; max = 0; mean = 0; stdDev = 0; sum = 0; sumsqr = 0; count = 0; surface = 0; area = 0; perimeter = 0;
	centroid.x = 0; centroid.y = 0;
	distScale = 1;
	areaScale = 1;
	distMenu.value = "meters";
	areaMenu.value = "square meters";
	cbRedraw();
	setDefaultWhenClose = true;
	}  # end of OnActivate 
 
# Called when tool is deactivated (usually when switching to another tool).
func OnDeactivate () {
	setDefaultWhenClose = false;
	cbClose();
	}  # end of OnDeactivate