crisp_filter.sml

# Created by Dan Glasser
# 18July03
# crisp filter
# compute the principle components
# sharpen the first component with a summary filter
# recombine using the sharpend pc1 and the other pc bands

numeric answer = PopupYesNo("Assuming there are 4 input bands, is this ok?", 1);
if (answer==0)
	Exit();

numeric numbands = 4;

raster band1, band2, band3, band4;

# assumption is 4 bands
GetInputRaster(band1);
GetInputRaster(band2);
GetInputRaster(band3);
GetInputRaster(band4);

string band1name$ = band1.$Info.Name;
string band2name$ = band2.$Info.Name;
string band3name$ = band3.$Info.Name;
string band4name$ = band4.$Info.Name;

class MATRIX forwardPca, inversePca, inMatrix, outMatrix;

# outMat = (forward * inMat)
outMatrix = CreateMatrix(numbands, 1);
# use to calculate the transformation
inMatrix = CreateMatrix(numbands, 1);
     
forwardPca = CreateMatrix(numbands, numbands);
inversePca = CreateMatrix(numbands, numbands);

# calculate principle components
PrincipalComponents(forwardPca, inversePca, band1, band2, band3, band4);

string filename$ = GetOutputFileName("pca", "Enter an output file name for the Principal Component and Filtered result rasters", ".rvc");
CreateProjectFile(filename$, "pca result");

# Create the rasters to hold principal components
# Output values will contain negative numbers so use signed type.
# Using "64-bit float"
Raster pc1, pc2, pc3, pc4;
CreateRaster(pc1, filename$, "pc1", "result of pca: pc1", NumLins(band1), NumCols(band1), "64-bit float");
CreateRaster(pc2, filename$, "pc2", "result of pca: pc2", NumLins(band2), NumCols(band2), "64-bit float");
CreateRaster(pc3, filename$, "pc3", "result of pca: pc3", NumLins(band3), NumCols(band3), "64-bit float");
CreateRaster(pc4, filename$, "pc4", "result of pca: pc4", NumLins(band4), NumCols(band4), "64-bit float");

# all bands have equal number of lins/cols
numeric r, c;
# save the principle component rasters
for r = 1 to NumLins(band1)
{
	for c = 1 to NumCols(band1)
	{
		# set inMatrix values to input bands values
		SetMatrixItem(inMatrix, 0, 0, band1[r,c]);
		SetMatrixItem(inMatrix, 1, 0, band2[r,c]);
		SetMatrixItem(inMatrix, 2, 0, band3[r,c]);
		SetMatrixItem(inMatrix, 3, 0, band4[r,c]);
		# out = forward * in
		MultiplyMatrix(outMatrix, forwardPca, inMatrix);
		pc1[r,c] = GetMatrixItem(outMatrix, 0, 0);
		pc2[r,c] = GetMatrixItem(outMatrix, 1, 0);
		pc3[r,c] = GetMatrixItem(outMatrix, 2, 0);
		pc4[r,c] = GetMatrixItem(outMatrix, 3, 0);
	}
}

# create a 3 x 3 summary filter
array numeric summaryfilter[3,3];

numeric i,j;
for i=1 to 3
{
	for j=1 to 3
	{
		summaryfilter[i,j] = -1;
	}
}
summaryfilter[2,2] = 10;

# 3x3 summary filter
# [ -1 -1 -1 ]
# [ -1 10 -1 ]
# [ -1 -1 -1 ]

Raster sharppc1;
CreateRaster(sharppc1, filename$, "sharppc1", "sharpendpc1", NumLins(band1), NumCols(band1), "64-bit float");

sharppc1 = FocalFilter(pc1, summaryfilter);

# create the rasters to hold the sharpened bands
Raster sharpenedBand1, sharpenedBand2, sharpenedBand3, sharpenedBand4;
CreateRaster(sharpenedBand1, filename$, "sharp_" + band1name$, "result of sharpening " + band1name$ + ": sharp_" + band1name$, NumLins(band1), NumCols(band1), "64-bit float");
CreateRaster(sharpenedBand2, filename$, "sharp_" + band2name$, "result of sharpening " + band2name$ + ": sharp_" + band2name$, NumLins(band2), NumCols(band2), "64-bit float");
CreateRaster(sharpenedBand3, filename$, "sharp_" + band3name$, "result of sharpening " + band3name$ + ": sharp_" + band3name$, NumLins(band3), NumCols(band3), "64-bit float");
CreateRaster(sharpenedBand4, filename$, "sharp_" + band4name$, "result of sharpening " + band4name$ + ": sharp_" + band4name$, NumLins(band4), NumCols(band4), "64-bit float");

# compute sharpened bands using inverse pca matrix and sharpend pc1
for r = 1 to NumLins(band1)
{
	for c = 1 to NumCols(band1)
	{
		# set inMatrix values to input bands values
		SetMatrixItem(inMatrix, 0, 0, sharppc1[r,c]);
		SetMatrixItem(inMatrix, 1, 0, pc2[r,c]);
		SetMatrixItem(inMatrix, 2, 0, pc3[r,c]);
		SetMatrixItem(inMatrix, 3, 0, pc4[r,c]);
		# out = forward * in
		MultiplyMatrix(outMatrix, inversePca, inMatrix);
		sharpenedBand1[r,c] = GetMatrixItem(outMatrix, 0, 0);
		sharpenedBand2[r,c] = GetMatrixItem(outMatrix, 1, 0);
		sharpenedBand3[r,c] = GetMatrixItem(outMatrix, 2, 0);
		sharpenedBand4[r,c] = GetMatrixItem(outMatrix, 3, 0);
	}
}

# copy any georef objects if available
CopySubobjects(band1, pc1, "georef");
CopySubobjects(band1, pc2, "georef");
CopySubobjects(band1, pc3, "georef");
CopySubobjects(band1, pc4, "georef");
CopySubobjects(band1, sharppc1, "georef");
CopySubobjects(band1, sharpenedBand1, "georef");
CopySubobjects(band1, sharpenedBand2, "georef");
CopySubobjects(band1, sharpenedBand3, "georef");
CopySubobjects(band1, sharpenedBand4, "georef");

# do clean up
CloseRaster(band1);
CloseRaster(band2);
CloseRaster(band3);
CloseRaster(band3);
CloseRaster(pc1);
CloseRaster(pc2);
CloseRaster(pc3);
CloseRaster(pc4);
CloseRaster(sharppc1);
CloseRaster(sharpenedBand1);
CloseRaster(sharpenedBand2);
CloseRaster(sharpenedBand3);
CloseRaster(sharpenedBand4);
DestroyMatrix(forwardPca);
DestroyMatrix(inversePca);
DestroyMatrix(inMatrix);
DestroyMatrix(outMatrix);