princomp.sml

# princomp.sml
# Computes principle component transformation matrix for a set of rasters.
clear(); # clear the console
class RVC_RASTER Rred, Rgreen, Rblue;
class RVC_RASTER R1, R2, R3;
class MATRIX forwardMat, inverseMat, inMat, outMat;
forwardMat = CreateMatrix(3,3); # create the transformation matrices
inverseMat = CreateMatrix(3,3);  # to send to princ. comp. function
inMat = CreateMatrix(3,1);  # outMat = (forward * inMat) 
outMat = CreateMatrix(3,1); # use to calculate the transformation
# get three input rasters to do principal components on
GetInputRaster(Rred);   # get red input raster 
if ( RastType(Rred)  <> "8-bit unsigned" ) {
   PopupMessage( "Only 8 bit unsigned types are supported." );
   CloseRaster(Rred);      # close it
   }
GetInputRaster(Rgreen);   # get green input raster
if ( RastType(Rgreen)  <> "8-bit unsigned" ) {
   PopupMessage( "Only 8 bit unsigned types are supported." );
   CloseRaster(Rgreen);      # close it
   }
GetInputRaster(Rblue);   # get blue input raster
if ( RastType(Rblue)  <> "8-bit unsigned" ) {
   PopupMessage( "Only 8 bit unsigned types are supported." );
   CloseRaster(Rblue);      # close it
   }
# create the rasters to hold principal components
# output values will contain negative numbers so 
# allways use signed type
GetOutputRaster(R1, NumLins(Rred), NumCols(Rred), "64-bit float");
GetOutputRaster(R2, NumLins(Rred), NumCols(Rred), "64-bit float");
GetOutputRaster(R3, NumLins(Rred), NumCols(Rred), "64-bit float");
# get the tranformation matrices
PrincipalComponents(forwardMat, inverseMat, Rred, Rgreen, Rblue );
# print the forward matrix
print( "Forward matrix" );
numeric r, c, a;
for r = 0 to 2 {
	for c = 0 to 2 {
		a = GetMatrixItem( forwardMat, r, c );
		printf( "%8.4f ", a );
		}
	printf( "%s\n\n", " " );
	}
# print the inverse matrix
print( "Inverse matrix" );
for r = 0 to 2 {
	for c = 0 to 2 {
		a = GetMatrixItem( inverseMat, r, c );
		printf( "%8.4f ", a );
		}
	printf( "%s\n\n", " " );
	}
# now compute the principal component values 
# for each point in raster
numeric numColumns = NumCols(Rred);
numeric numLines = NumLins(Rred);
for r = 1 to NumLins(Rred) {
	for c = 1 to NumCols(Rred) {
		SetMatrixItem( inMat, 0, 0, Rred[r,c] );
		SetMatrixItem( inMat, 1, 0, Rgreen[r,c] );
		SetMatrixItem( inMat, 2, 0, Rblue[r,c] ) ;
		MultiplyMatrix( outMat, forwardMat, inMat );
		R1[r,c] = GetMatrixItem( outMat, 0, 0 );
		R2[r,c] = GetMatrixItem( outMat, 1, 0 );
		R3[r,c] = GetMatrixItem( outMat, 2, 0 );
		}
	}
# do clean up
CloseRaster(Rred);
CloseRaster(Rgreen);
CloseRaster(Rblue);
CloseRaster(R1);
CloseRaster(R2);
CloseRaster(R3);
DestroyMatrix( forwardMat );
DestroyMatrix( inverseMat );
DestroyMatrix( inMat );
DestroyMatrix( outMat );