Generate parameters for a world file without using the REGISTER command

/*
/* AFFINE.AML -- Compute Affine Transformation Parameters
/*
/* Macro for computing the parameters that describe the transformation from
/* one set of coordinate points (the source coordinates) to another set (the
/* destination coordinates). Six parameters are calculated, and these are
/* described in the "ARC/INFO Image Integrator User's Guide", Chapter 2. The
/* parameters are then used to calculate transformed destination coordinates,
/* from which error residuals and a total destination root mean square (RMS)
/* error value are calculated.
/*
/* Input: .TRN$OSX original source coordinates
/* .TRN$OSY
/* .TRN$ODX original destination coordinates
/* .TRN$ODY
/*
/* Output: .TRN$XP1 X transformation parameter 1 (C in the User's Guide)
/* .TRN$XP2 X transformation parameter 2 (A)
/* .TRN$XP3 X transformation parameter 3 (B)
/* .TRN$YP1 Y transformation parameter 1 (F)
/* .TRN$YP2 Y transformation parameter 2 (D)
/* .TRN$YP3 Y transformation parameter 3 (E)
/* .TRN$TDX transformed destination coordinates
/* .TRN$TDY
/* .TRN$RDX error residuals (transformed - original)
/* .TRN$RDY
/* .TRN$DRMS destination root mean square error
/*
/*============================================================================

&severity &error &routine error
&args in_x in_y out_x out_y

/* See if you need to display usages.

&if [null %out_y%] &then &goto usage
&if [null %out_x%] &then &goto usage
&if [null %in_y%] &then &goto usage
&if [null %in_x%] &then &goto usage


/* -- Initialize the matrices; TOP and BOT are 1x3 arrays, and TMATRIX is
/* -- a 3x3 array.

&do i = 1 &to 3
&s top%i% = 0
&s bot%i% = 0
&do j = 1 &to 3
&s tmatrix%i%%j% = 0
&end
&end

/* -- Load the matrices with the coordinates.

&s i = 1
&do &while [variable .trn$osx%i%]

/* -- Ensure all coordinate variables are present for this pass through
/* -- the loop.

/*&if not [variable .trn$osy%i%] or ~
/* not [variable .trn$odx%i%] or ~
/* not [variable .trn$ody%i%] &then &do
/* &type ERROR: Missing original transformation coordinate variables (.TRN$O*)
/* &call error
/* &end

/* -- AML loses precision when dividing by very large numbers; subtract
/* -- the first destination coordinates from all subsequent destination
/* -- coordinates, translating those points so that the first point lies on
/* -- the origin of the Cartesian plane; the first destination coordinates
/* -- will be added back into the transformation parameters later.

&s odx%i% = [value .trn$odx%i%] - %.trn$odx1%
&s ody%i% = [value .trn$ody%i%] - %.trn$ody1%

/* -- Add the coordinates to the transformation matrix TMATRIX.

&s tmatrix11 = %tmatrix11% + 1
&s tmatrix12 = %tmatrix12% + [value .trn$osx%i%]
&s tmatrix13 = %tmatrix13% + [value .trn$osy%i%]
&s tmatrix22 = %tmatrix22% + ( [value .trn$osx%i%] * [value .trn$osx%i%] )
&s tmatrix23 = %tmatrix23% + ( [value .trn$osx%i%] * [value .trn$osy%i%] )
&s tmatrix33 = %tmatrix33% + ( [value .trn$osy%i%] * [value .trn$osy%i%] )

/* -- Add the coordinates to the TOP and BOT matrices.

&s top1 = %top1% + [value odx%i%]
&s top2 = %top2% + ( [value .trn$osx%i%] * [value odx%i%] )
&s top3 = %top3% + ( [value .trn$osy%i%] * [value odx%i%] )
&s bot1 = %bot1% + [value ody%i%]
&s bot2 = %bot2% + ( [value .trn$osx%i%] * [value ody%i%] )
&s bot3 = %bot3% + ( [value .trn$osy%i%] * [value ody%i%] )

&s i = %i% + 1
&end

&s tmatrix21 = %tmatrix12%
&s tmatrix32 = %tmatrix23%
&s tmatrix31 = %tmatrix13%

/* -- Perform matrix algebra; sorry, but this is a "black box".
/* -- The code is from a program in
Visual Basic.

&do i = 1 &to 3
&do j = 1 &to 3
&if %j% ne %i% &then &do
&s tmatrix%i%%j% = [value tmatrix%i%%j%] / [value tmatrix%i%%i%]
&end
&end
&s tmatrix%i%%i% = 1 / [value tmatrix%i%%i%]
&do j = 1 &to 3
&if %j% ne %i% &then &do
&do k = 1 &to 3
&if %k% ne %i% &then &do
&s tmatrix%j%%k% = [value tmatrix%j%%k%] - ~
( [value tmatrix%j%%i%] * [value tmatrix%i%%k%] )
&end
&end
&s tmatrix%j%%i% = ( 0 - [value tmatrix%j%%i%] ) * [value tmatrix%i%%i%]
&end
&end
&end

/* -- Calculate transformation parameters from the matrix.

&do i = 1 &to 3
&s .trn$xp%i% = 0
&s .trn$yp%i% = 0
&do j = 1 &to 3
&s .trn$xp%i% = [value .trn$xp%i%] + ~
( [value tmatrix%i%%j%] * [value top%j%] )
&s .trn$yp%i% = [value .trn$yp%i%] + ~
( [value tmatrix%i%%j%] * [value bot%j%] )
&end
&end

/* -- The first X and Y parameters describe the translation for the
/* -- transformation; add the first destination coordinates, which were
/* -- subtracted prior to the matrix algebra, to these parameters.

&s .trn$xp1 = %.trn$xp1% + %.trn$odx1%
&s .trn$yp1 = %.trn$yp1% + %.trn$ody1%

/* -- Using the transformation parameters, calculate adjusted destination
/* -- coordinates and error residuals.

&s .trn$drms = 0
&s i = 1
&do &while [variable .trn$odx%i%]
&s .trn$tdx%i% = %.trn$xp1% + ~
( [value .trn$osx%i%] * %.trn$xp2% ) + ~
( [value .trn$osy%i%] * %.trn$xp3% )
&s .trn$tdy%i% = %.trn$yp1% + ~
( [value .trn$osx%i%] * %.trn$yp2% ) + ~
( [value .trn$osy%i%] * %.trn$yp3% )
&s .trn$rdx%i% = [value .trn$tdx%i%] - [value .trn$odx%i%]
&s .trn$rdy%i% = [value .trn$tdy%i%] - [value .trn$ody%i%]

/* -- Accumulate the squared distance between the original and transformed
/* -- destination points.

&s dist = [invdistance 0 0 [value .trn$rdx%i%] [value .trn$rdy%i%]]
&s .trn$drms = %.trn$drms% + ( %dist% * %dist% )

&s i = %i% + 1
&end
&s i = %i% - 1

/* -- Calculate the root mean square error.

&s .trn$drms = %.trn$drms% / ( %i% - 1 )
&s .trn$drms = [sqrt %.trn$drms%]

&return

/* Display usage.

&label usage
&type Usage: AFFINE <out_y)
&return

/*===========================================================================

&routine error
&severity &error &ignore
&return; &return &error Bailing out of AFFINE.AML ...