clear();
#
# reclassify.sml
#
# This script is written to reclassify cell values in the user-entered input raster.
# Reclassification is based on meeting certain criteria for class (value) and area.
# For modification purposes, the variables used for this purpose are "threshold",
# "oldValue", "and newValue".
# Currently it is set to reclassify all areas of <5 acres in class '128' to class '0'.
#

class RASTER Rin;
class RASTER Result;
class STATUSDIALOG statusdialog;
class STATUSCONTEXT statuscontext;

# create the 5 acre threshold
numeric thresholdAcres = 5;

numeric oldValue = 128, newValue = 0;

# Get the input and output rasters.
GetInputRaster(Rin);
IgnoreNull(Rin);
GetOutputRaster(Result, NumLins(Rin), NumCols(Rin), RastType(Rin));
IgnoreNull(Result);

# Copy all cell values
printf("Copying source raster\n");
Result = Rin;

# Get the georeference info
CopySubobjects(Rin, Result);	# get colormap, contrast, georef, etc.

# Convert the threshold from Acres to cells
numeric threshold = thresholdAcres * 4046.85642 / (LinScale(Rin) * ColScale(Rin));

printf("Eliminating solid features less than or equal to %g acres (%d cells)\n", thresholdAcres, threshold);

numeric lin, col, r, up, left, count;

foreach Result[lin, col] {
	r = Result[lin,col];
	if (r == oldValue) {
		up = -1;
		left = -1;
		if (col > 1) left = Result[lin, col-1];
		if (lin > 1) up = Result[lin-1, col];

		#printf("R[%d,%d] = %d,  Cell to the left = %d, Cell above = %d\n", lin, col, r, left, up);

		if (left != oldValue && up != oldValue) {
			#printf("Checking at %d,%d\n", lin, col);
			count = RasterGetSolidAreaSize(Result, lin, col, threshold+1);
			if (count <= threshold) {
				#printf("Filling at %d,%d - count = %d\n", lin, col, count);
				RasterFloodFill(Result, lin, col, newValue);
				}
			}
		}
	}

CreateHistogram(Result);
CreatePyramid(Result);