biomass.sml

#
#	This is how to declare a global variable so that
#	they can be used in functions before it's assigned.
#	Without the explicit declaration, the variable will
#	be assumed to be numeric or a Raster (depending on
#	the case of the first letter).
#
class RVC_COLORMAP cmap;
class COLOR black, white, c0;
class XmDrawingArea da;
class XmForm form;			# The main form dialog
class GC gc;	# The context for drawing into the view
raster Composite, BiomassRaster, FilteredBiomassRaster;
vector BiomassVector;
class GRE_VIEW view;
class GRE_GROUP group;
class GRE_LAYER_RASTER layer, overlay;			# Our base layer
class GRE_LAYER_VECTOR voverlay;
class GRE_TOOL_POLYLINE tool;
class tooltip tip;
region reg;
array numeric histo[256];
numeric cellsize, unitconv;
numeric doingtest;
numeric haveBiomassRaster, haveBiomassVector;
class XmToggleButton drawhisto;
class PushButtonItem buttonConvertToVector, buttonFilterVector;
class PushButtonItem buttonClearOverlay, buttonExit;
class XmForm button_row;
class XmFrame frame;
numeric height, width;

unitconv = GetUnitConvArea("square meters", "acres");
numeric version = round(_context.Version * 10);
if (version < 74) {
	PopupMessage("This version requires that you are using TNTmips2008:74 or later");
	Exit();
	}


proc setTip () {
	if (!haveBiomassRaster) {
		if (tool.HasPosition) {
			tip.string = "Press right button\nto close polygon";
			}
		else {
			tip.string = "Press the left button\nand draw around a field.";
			}
		}
	else {
		tip.string = "";
		}
	}


#
#	This is the callback function which will be called
#	when button1 is pressed.  All widget callbacks get
#	two parameters (for now).  The first is the widget
#	generating the callback.  The second is an optional
#	"cbdata" and will be whatever was passed to  
#	WidgetAddCallback for the last parameter.  It isn't
#	really supported yet.
#	When done, callbacks will also get a 3rd parameter
#	which is the Callback Struct that the XmWidget passes
#	There's currently no way to declare these things
#

proc cbQuit(class widget widget) {
	DialogClose(form);
	}

numeric sample_height, sample_width, sample_space, sample_text_x, num_ranges;

sample_height = 15;
sample_width = 40;
sample_space = 5;
sample_text_x = 37;
num_ranges = 10;
string str$ = "";

#
#	Draw the bar graph
#
proc drawGraph(class XmDrawingArea da) {
	local numeric maxval, total, i, scale, w, y;

	# Clear the drawing area
	ActivateGC(gc);
	SetColorRGB(0,0,0);
	FillRect(0, 0, da.width, da.height);

	maxval = 0;
	total = 0;
	for i = 1 to num_ranges {
		maxval = SetMax(maxval, histo[i]);
		total += histo[i];
		}
	if (total) {	# Don't draw bars if nothing selected
		w = da.width - 2 * sample_space;
		scale = w / maxval;
		if (drawhisto.set) {
			for i = 1 to num_ranges {
				y = (num_ranges - i + 1) * (sample_height + sample_space) + sample_space;
				SetColor(ColorMapGetColor(cmap, i));
				w = scale * histo[i];
				if (w < 1 && histo[i] != 0) w = 1;
				FillRect(sample_space, y, w, sample_height);
				}
			}
		else {
			DrawTextSetColors(white, black);
			y = sample_space;
			gc.TextStyle.JustifyRight = 0;
			gc.TextStyle.JustifyCenter = 1;
			DrawTextSimple("Acres", da.width / 2, y + sample_height - 4);
			gc.TextStyle.JustifyRight = 1;
			gc.TextStyle.JustifyCenter = 0;
			DrawTextSetColors(black, white);
			for i = 1 to num_ranges {
				y = (num_ranges - i + 1) * (sample_height + sample_space) + sample_space;
				SetColor(ColorMapGetColor(cmap, i));
				FillRect(sample_space, y, sample_width, sample_height);
				SetColorName("white");
				DrawRect(sample_space, y, sample_width, sample_height);
				str$ = sprintf("%.2f", histo[i] * cellsize * unitconv);
				DrawTextSimple(str$, sample_text_x, y + sample_height - 4);
				}
			y = (num_ranges + 1) * (sample_height + sample_space) + sample_space;
			str$ = sprintf("%.2f", total * cellsize * unitconv);
			DrawTextSimple(str$, sample_text_x, y + sample_height - 4);
			}
		}
	}


#
#	Callback for when user clicks the right mouse button
#	on the polygon tool (or hits apply)
#
proc cbToolApply(class GUI_GADGET_REGION tool) {
	local numeric i, val, scale, minval, maxval, ir, red, total, offset, lo, hi;

	reg = tool.Region;
	FillRegion(RegionTrans(reg, ViewGetTransViewToScreen(view)));
	reg = RegionTrans(reg, ViewGetTransLayerToView(view, layer, 1));
	StatusSetDefaultHandle(view.StatusHandle);

	ViewSetMessage(view, "Filling raster with 0");
	BiomassRaster = 0;
	if (1) {
		ViewSetMessage(view, "Computing scale");
		for i=0 to 200 histo[i] = 0;
        for each Composite in reg.$Data {
			ir = Composite.red;
			red = Composite.green;
			val = ((ir - red) / (ir + red) + 1) * 100;
			histo[val] = histo[val] + 1;
			}
		total = 0;
		for i=0 to 200 total += histo[i];
		lo = total * .02;
		hi = total * .98;
		minval = -100;
		maxval = -100;
		total = 0;
		for i = 0 to 200 {
			if (total > lo && minval < -1) {
				minval = (i / 100) - 1;
				}
			total += histo[i];
			if (total > hi && maxval < -1) {
				maxval = (i / 100) - 1;
				}
			}
        scale = num_ranges / (maxval - minval);
        offset = -minval;
        printf("offset = %f, scale = %f\n", offset, scale);
		}
	else {
		scale = num_ranges / 2.0;
		offset =  1.0;
		}

	ViewSetMessage(view, "Processing region");
	for i = 1 to num_ranges histo[i] = 0;
	for each Composite in reg.$Data {
		ir = Composite.red;
		red = Composite.green;
		val = (ir - red) / (ir + red);
		val = Bound((val + offset) * scale + 1, 1, num_ranges);
		histo[val] = histo[val] + 1;
		BiomassRaster = val;
		}
	CloseRaster(BiomassRaster);
	LayerShow(overlay, view);
	ViewRedraw(view);
	drawGraph(da);

	# Do something with the region here
	haveBiomassRaster = 1;
	tool.HasPosition = 0;	# Clears the current position
	buttonConvertToVector.Disabled = 0;
	setTip();
	}

proc cbConvToVector () {
	StatusSetDefaultHandle(view.StatusHandle);

	if (0) {		# Make this 1 to low-pass filter the raster first
		FilteredBiomassRaster = 0;
		ViewSetMessage(view, "Lowpass filter");
		foreach BiomassRaster in reg.$Data {
			FilteredBiomassRaster = FocalMedian(BiomassRaster, 3, 3);
			}
		ViewSetMessage(view, "Converting Raster to Vector");
		BiomassVector = RasterToVectorBound(FilteredBiomassRaster);
		}
	else {
		ViewSetMessage(view, "Converting Raster to Vector");
		BiomassVector = RasterToVectorBound(BiomassRaster);
		}

	CloseVector(BiomassVector);
	if (!voverlay) {
		voverlay = GroupQuickAddVectorVar(group, BiomassVector);
		voverlay.Poly.NormalStyle.BorderColor.name = "black";
		voverlay.Line.NormalStyle.Color.name = "black";
		}
	LayerShow(voverlay, view);
	ViewRedraw(view);
	buttonFilterVector.Disabled = 0;
	}

proc cbFilterVector () {
	local numeric i, numPolys, numToDelete;

	FilteredBiomassRaster = 0;
	StatusSetDefaultHandle(view.StatusHandle);
	ViewSetMessage(view, "Generating std stats table");
	VectorToolkitInit(BiomassVector);
	VectorUpdateStdAttributes(BiomassVector);		# Compute std stats table
	numToDelete = 0;
	numPolys = NumVectorPolys(BiomassVector);
	array numeric deleteList[numPolys];
	ViewSetMessage(view, "Looking for polygons to delete");
	for i=1 to numPolys {
		if (BiomassVector.poly[i].POLYSTATS.Area < 40) {
			numToDelete += 1;
			deleteList[numToDelete] = i;
			}
		}
	if (numToDelete) {
		printf("Found %d polygons to delete\n", numToDelete);
		ViewSetMessage(view, "Deleting tiny polygons");
		# Now delete std attributes do delete polys doesn't have to maintain them
		VectorSetFlags(BiomassVector, "NoDBStatTable");
		VectorDeleteStdAttributes(BiomassVector);
#		VectorDeletePolys(BiomassVector, deleteList, numToDelete, "LongestLine");
		VectorDeletePolys(BiomassVector, deleteList, numToDelete);
		}

	ViewSetMessage(view, "Closing vector");
	CloseVector(BiomassVector);
	ViewRedraw(view);
	}

proc cbHideOverlay(class widget widget) {
	local numeric i;
	LayerHide(overlay, view);
	for i = 0 to 10 histo[i] = 0;
	ViewRedraw(view);
	drawGraph(da);
	}

proc cbRedrawGraph() {
	drawGraph(da);
	}

proc cbToolPosSet(class GUI_GADGET tool) {
	setTip();
	}

#####################################################
#
#	MAIN PROGRAM
#

# Get an input raster.  We will use it's size to determine 
# 	 how big to make our window
GetInputRaster(Composite);


#
#	Create the Biomass Raster that we'll use for the overlay
#
string outfilename$ = "biomass.rvc";
DeleteFile(outfilename$);
CreateRaster(BiomassRaster, outfilename$, "biomass", "Biomass output", NumLins(Composite), NumCols(Composite), "8-bit unsigned");
#CreateTempRaster(BiomassRaster, NumLins(Composite), NumCols(Composite), "8-bit unsigned");
#GetOutputRaster(BiomassRaster, NumLins(Composite), NumCols(Composite), "8-bit unsigned");
CopySubobjects(Composite, BiomassRaster, "georef");
cellsize = LinScale(Composite) * ColScale(Composite); 
c0 = ColorMapGetColor(cmap, 0);
c0.transp = 50;		# Make it 50% transparent
ColorMapSetColorHIS(cmap, 10, 230, 50, 100);	# Green
ColorMapSetColorHIS(cmap, 9, 180, 70, 100);		# Yellow
ColorMapSetColorHIS(cmap, 8, 150, 70, 100);		# Lt Orange
ColorMapSetColorHIS(cmap, 7, 144, 50, 100);		# Dk Orange
ColorMapSetColorHIS(cmap, 6, 120, 50, 100);		# Red
#ColorMapSetColorHIS(cmap, 5, 96, 50, 100);
ColorMapSetColorHIS(cmap, 5, 72, 50, 100);		# Magenta
ColorMapSetColorHIS(cmap, 4, 48, 50, 100);		# Purple
ColorMapSetColorHIS(cmap, 3, 10, 50, 100);		# Blue
ColorMapSetColorHIS(cmap, 2, 342, 65, 100);		# Lt Blue
ColorMapSetColorHIS(cmap, 1, 320, 70, 100);		# Cyan
ColorMapSetColor(cmap, 0, c0);
ColorMapWriteToRastVar(BiomassRaster, cmap, "ColorMap");
SetNull(BiomassRaster, 0);


CreateTempRaster(FilteredBiomassRaster, NumLins(Composite), NumCols(Composite), "8-bit unsigned");
CopySubobjects(Composite, FilteredBiomassRaster, "georef");


#
#	Create a temp vector for the biomass->vector conversion
#
#CreateTempVector(BiomassVector);
CreateVector(BiomassVector, outfilename$, "BiomassVect", "Biomass vector");


#  Create a modal dialog (aka "window").  A modal dialog
#   prevents other windows from accepting input until it
#   goes away.  This is the easiest kind to work with.
#	Non-modal dialogs will be supported later
#
#	This function returns a form that can be used as a
#	parent for other widgets
#
#	If we had another dialog which we wanted to use as
#	a "parent" to this one, we would pass it as the
#	2nd parameter to CreateModalFormDialog().

form = CreateModalFormDialog("Relative Biomass Mapping");

#
#	Create a push button on the form.  Set its attachments
#	so that sets in the lower left of the form.  To attach
#	an edge to the form rather than another widget, just
#	set it to the form.  In Motif, this would be done by
#	setting the attachment to XmATTACH_FORM, which I plan
#	to implement here too.  Setting the RightPosition to 50
#	ties the right edge of the button to a spot 50% of the
#	width of the form, (assuming form.FractionBase still
#	has the default value of 100) no matter how wide the 
#	form gets
#

buttonClearOverlay = CreatePushButtonItem("Clear Overlay", cbHideOverlay);
buttonConvertToVector = CreatePushButtonItem("Convert To Vector", cbConvToVector);
buttonFilterVector = CreatePushButtonItem("Filter Vector", cbFilterVector);
buttonExit = CreatePushButtonItem("Exit", cbQuit);
buttonConvertToVector.Disabled = 1;
buttonFilterVector.Disabled = 1;
button_row = CreateButtonRow(form, buttonClearOverlay, buttonConvertToVector, buttonFilterVector, buttonExit);
button_row.BottomWidget = form;

drawhisto = CreateIconToggleButton(form, "RVCobjects", "histo", "Show as histogram");
drawhisto.LeftWidget = form;
drawhisto.TopWidget = form;
WidgetAddCallback(drawhisto.ValueChangedCallback, cbRedrawGraph);

frame = CreateFrame(form);
frame.TopWidget = drawhisto;
frame.BottomWidget = button_row;
frame.LeftWidget = form;

da = CreateDrawingArea(frame, 180, sample_width + 2 * sample_space + 2);
WidgetAddCallback(da.ExposeCallback, cbRedrawGraph);

#
#	Create a 2D group to display the raster in and quick-add it.
#
group = DispCreate2DGroup();
layer = GroupQuickAddRasterVar(group, Composite);
overlay = GroupQuickAddRasterVar(group, BiomassRaster);
overlay.NullCellsTransparent = 0;
overlay.UsesTransparency = 1;
overlay.DoBlendMask = 1;

#
#	Create a view for this group in our form.  Default its size to
#	the size of the raster.  For a real program, you'd want to
#	limit this to something reasonable.  
#

height = NumLins(Composite);
width = NumCols(Composite);
while (height > 400 or width > 640) {
	height /= 2;
	width /= 2;
	}

view = GroupCreateView(group, form, "view", height, width);
#
#	Make widget attachments for the view.  Connect it to the form
#	on all edges except the bottom.  There, attach it to one of the
#	buttons.  
#	Some notes about making attachment:
#		1) If you attach widget A to Widget B, don't attach B to A
#		   or anything else which would cause a circular reference.
#		2) The widget created and attached last is what will resize
#		   with the window.  If I had created the view first and then
#		   the buttons, resizing the window would make the buttons
#		   bigger.


form.BottomWidget = button_row;


LayerHide(overlay, view);
view.ShowDataTips = 0;

# Hide the whole icon bar to show that we can do it.
# (commented out because I don't really want it hidden)
#view.IconBar.hidden = 1;

#
#	Hide the scale position.  Do this BEFORE adding the
#	EndDrawCallback below. Turning off the ScalePos row
#	Will cause a redraw.  If the callback function tries
#	to draw onto a drawing area that doesn't exist yet,
#	it will IllOp.  (Yes, now you too can program your
#	own IllOps!)
#
#view.ScalePosVisible = 0;


#
#	Add a callback to the view's DrawEndCallback.
#	This will be called when the drawing is done.
#
#	Before I'm done, I plan to merge WidgetAddCallback()
#	and DispAddCallback() into a single AddCallback().
#	This will be done by making all callback lists a
#	subclass of a single class.  
#
DispAddCallback(view.RestoreCallback, cbRedrawGraph);

#
#	Create a region tool on the view
#
tool = ViewCreatePolygonTool(view, "My Tool", "edit_query", "standard");
ToolAddCallback(tool.ActivateCallback, cbToolApply);
ToolAddCallback(tool.PositionSetCallback, cbToolPosSet);
ToolAddCallback(tool.PositionClearedCallback, cbToolPosSet);

#
#	Have to call this after adding any tools of our own
#
#ViewAddToolIcons(view);
ViewAddStandardTools(view);

#
#	Open the dialog Need to do this before creating the GC
#	because until then the drawing area has no window.
#	Would like a way around this.  Could my CreateGC
#	realize it first if necessary?
#
DialogOpen(form);

#
#	Create a GC (Graphics Context) for the drawing area
#	You can have multiple gc's per drawing area.  Drawing
#	functions will always use the active GC.  
#
gc = CreateGCForDrawingArea(da);
gc.TextStyle.JustifyRight = 1;
# Clear the drawing area
ActivateGC(gc);
SetColorRGB(0,0,0);
black.name = "black";
white.name = "white";
FillRect(0, 0, da.width, da.height);
DrawTextSetFont("stork.of");
DrawTextSetColors(black, white);
DrawTextSetHeightPixels(9);

ViewRedraw(view);
#ViewActivateTool(view, tool);
tool.Managed = 1;
ViewSetMessage(view, "Draw around a field and press \"Apply\"");
tip = CreateToolTip(view.DrawingArea, "Draw around a field\nand press \"Apply\"");
tip.Disabled = 0;
tip.Delay = 5000;	# 5 seconds (5000 ms)

setTip();

#StatusSetDefaultHandle(view.StatusHandle);
#
#	Wait for the dialog to close.  This function will
#	not return until the dialog closes.  
#	The normal way to do this is have a button or menu
#	item who's callback calls DialogClose(form).
#	Clicking on the "X" button on the title bar will 
#	close it (although I plan to let you override that 
#	too)
#
DialogWaitForClose(form);

#
#	Destroy our form to free up its memory.  Destroying the
#	main form will also destroy all its children.
#	Need to destroy the view before the form because
#	otherwise the tool gets destroyed before the view
#	and the view tries to disable a nonexistant tool.
#	Result: IllOP!
#
ViewDestroy(view);
GroupDestroy(group);
DestroyWidget(form);