# Created by Dan Glasser
# 18July03
# Computes principal components of 4 input bands

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

numeric numbands = 4;

class RVC_RASTER band1, band2, band3, band4;

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

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 result", ".rvc");
CreateProjectFile(filename$, "pca result");

# create the rasters to hold principal components
# output values will contain negative numbers so
# always use signed type
class RVC_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;
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);
	}
}

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