# REGSTATS.SML - Allows user to select an area of a non-composite raster.
# 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 6.4

# 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;
class XmSeparator line1, line2;
class MdispRegionTool tool;
class GRE_LAYER rasterLayer;
class Raster targetRaster;
class XmDrawingArea da;
class GC gc;
string rasterName$;
class PushButtonItem saveButton, closeButton;
class XmOptionMenu distMenu, areaMenu;
numeric min, max, mean, stdDev, count, cells, sum, sumsqr, surface, area, perimeter;
class POINT2D centroid;
numeric distScale, areaScale;
class GRE_GROUP group;
numeric setDefaultWhenClose = false;

# Checks layer to see if it is valid.
func checkLayer() {
	local numeric valid = false;
	# Get name of active layer if it is usable.  If not output an error message.
	if (group.ActiveLayer.TypeID == "") {
		rasterName$ = "Group has no layers!";
		}
	else if (group.ActiveLayer.TypeID == "Raster") {
		rasterLayer = group.ActiveLayer;
		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!";
			}
		else {
			rasterName$ = rasterLayer.Name;
			valid = true;
			}
		}
	else
		rasterName$ = "Not a raster!";
	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");
	DrawTextSetFont("arial.ttf");
	DrawTextSetHeight(14);
	if (cells > 0) {
		DrawTextSimple(sprintf("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), 0, 10);
		}
	else DrawTextSimple(sprintf("Raster: %s\nCells:\nNull Cells:\nMinimum:\nMaximum:\nMean:\nStandard Deviation:\nArea:\nPerimeter:\nCentroid:\nSurface Area:", 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 RegionTool 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;
		local class REGION2D MyRgn;
		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.ActiveLayer.zscale;

		# Get region and transform it to the appropriate coordinate system.
		MyRgn = tool.RegionData;
		MyRgn = RegionTrans(MyRgn, ViewGetTransViewToScreen(View, 1));
		MyRgn = RegionTrans(MyRgn, ViewGetTransMapToView(View, rasterLayer.Projection, 1));

		# 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.
		numeric lin, col;
		foreach targetRaster[lin, col] in MyRgn {
			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 = MyRgn.GetArea();
			perimeter = MyRgn.GetPerimeter();
			centroid = MyRgn.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 clse button is pressed.  Closes the dialogs.
proc cbClose() {
	tool.Managed = 0;
	DialogClose(form);
	if (setDefaultWhenClose) {
		setDefaultWhenClose = false;
		View.SetDefaultTool();
		}
	}

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

	# Drawing area for displaying text.
	da = CreateDrawingArea(form, 173, 301);
	da.topWidget = form;
	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;
	buttonRow.bottomWidget = form;

	# Add create polygon tool.
	tool = ViewCreatePolygonTool(View);
	ToolAddCallback(tool.ActivateCallback, cbToolApply);
	}  # end of OnInitialize 
 
# Called when tool is to be destroyed, will not be called if tool was never activated.
func OnDestroy () {
	tool.Managed = 0;
	DestroyGC(gc);
	DestroyWidget(form);
	}  # end of OnDestroy 
 
# Called when tool is activated.
func OnActivate () {
	checkLayer();
	tool.Managed = 1;
	tool.HasPosition = 0;
	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