path: root/pkg/images/immatch/doc/geomap.hlp

.help geomap Jan01 images.immatch
.ih
NAME
geomap -- compute one or more spatial transformation functions
.ih
USAGE
geomap input database xmin xmax ymin ymax
.ih
PARAMETERS
.ls input
The list of text files containing the pixel coordinates of control points in
the reference and input images. The control points are listed
one per line with xref, yref, xin, and yin in columns 1 through 4 respectively.
.le
.ls database
The name of the text file database where the computed transformations will
be stored.
.le
.ls xmin, xmax, ymin, ymax
The range of reference coordinates over which the computed coordinate
transformation is valid. If the user is working in pixel units  these limits
should normally be set to the values of the column and row limits of the
reference image, e.g xmin = 1.0, xmax = 512, ymin= 1.0, ymax = 512 for
a 512 x 512 image. The minimum and maximum xref and yref values in \fIinput\fR
are used if xmin, xmax, ymin, or ymax are undefined.
.le
.ls transforms = ""
An optional list of transform record names. If transforms is undefined 
the database record(s) are assigned the names of the
individual text files specified by \fIinput\fR.
.le
.ls results = ""
Optional output files containing a summary of the results including a
description of the transform geometry and a listing of the input coordinates,
the fitted coordinates, and the fit residuals. The number of results files
must be one or equal to the number of input files. If results is "STDOUT" the
 results summary is printed on the standard output.
.le
.ls fitgeometry = "general"
The fitting geometry to be used. The options are the following.
.ls shift
X and y shifts only are fit.
.le
.ls xyscale
X and y shifts and x and y magnification factors are fit. Axis flips are
allowed for.
.le
.ls rotate
X and y shifts and a rotation angle are fit. Axis flips are allowed for.
.le
.ls rscale
X and y shifts, a magnification factor assumed to be the same in x and y, and a
rotation angle are fit. Axis flips are allowed for.
.le
.ls rxyscale
X and y shifts, x and y magnifications factors, and a rotation angle are fit.
Axis flips are allowed for.
.le
.ls general
A polynomial of arbitrary order in x and y is fit. A linear term and a
distortion term are computed separately. The linear term includes an x and y
shift, an x and y scale factor, a rotation and a skew.  Axis flips are also
allowed for in the linear portion of the fit. The distortion term consists
of a polynomial fit to the residuals of the linear term. By default the
distortion term is set to zero.
.le

For all the fitting geometries except "general" no distortion term is fit,
i.e. the x and y polynomial orders are assumed to be 2 and the cross term
switches are assumed to be "none", regardless of the values of the
\fIxxorder\fR, \fIxyorder\fR, \fIxxterms\fR, \fIyxorder\fR, \fIyyorder\fR and
\fIyxterms\fR parameters set by the user.
.le
.ls function = "polynomial"
The type of analytic surface to be fit. The options are the following.
.ls legendre
Legendre polynomials in x and y.
.le
.ls chebyshev
Chebyshev polynomials in x and y.
.le
.ls polynomial
Power series in x and y.
.le
.le
.ls xxorder = 2, xyorder = 2,  yxorder = 2, yyorder = 2
The order of the polynomials in x and y for the x and y fits respectively.
The default order and cross term settings define the linear term in x
and y, where the 6 coefficients can be interpreted in terms of an x and y shift,
an x and y scale change, and rotations of the x and y axes. The "shift",
"xyscale", "rotation", "rscale", and "rxyscale", fitting geometries
assume that the polynomial order parameters are 2 regardless of the values
set by the user. If any of the order parameters are higher than 2 and
\fIfitgeometry\fR is "general", then a distortion surface is fit to the
residuals from the linear portion of the fit.
.le
.ls xxterms = "half", yxterms = "half"
The options are:
.ls none
The individual polynomial terms contain powers of x or powers of y but not
powers of both.
.le
.ls half
The individual polynomial terms contain powers of x and powers of y, whose
maximum combined power is max (xxorder - 1, xyorder - 1) for the x fit and
max (yxorder - 1, yyorder - 1) for the y fit. 
.le
.ls full
The individual polynomial terms contain powers of x and powers of y, whose
maximum combined power is max (xxorder - 1, xyorder - 1) for the x fit and
max (yxorder - 1, yyorder - 1) for the y fit.
.le

The "shift", "xyscale", "rotation", "rscale", and "rxyscale" fitting
geometries, assume that the cross term switches are set to "none"
regardless of the values set by the user.  If either of the cross terms
parameters are set to "half" or "full" and \fIfitgeometry\fR is "general"
then a distortion surface is fit to the residuals from the linear
portion of the fit.
.le
.ls maxiter = 0
The maximum number of rejection iterations. The default is no rejection.
.le
.ls reject = 3.0
The rejection limit in units of sigma.
.le
.ls calctype = "real"
The precision of the coordinate transformation calculations. The options are
real and double.
.le
.ls verbose = yes
Print messages about actions taken by the task ?
.le
.ls interactive = yes
In interactive mode the user may interact with the fitting process, e.g.
change the order of the fit, reject points, display the data, etc.
.le
.ls graphics = "stdgraph"
The graphics device.
.le
.ls cursor = ""
The graphics cursor.
.le
.ih
DESCRIPTION

GEOMAP computes the transformation required to map the reference coordinate
system to the input coordinate system.  The coordinates of points in common
to the two systems are listed in the input text file(s) \fIinput\fR
one per line in the following format: "xref yref xin yin".

The computed transforms are stored in the text database file \fIdatabase\fR
in records with names specified by the parameter \fItransforms\fR. If the
transforms parameter is undefined the records are assigned the name of
the input coordinate files.

The computed transformation has the form shown below, where the reference
coordinates must be defined in the coordinate system of the reference image
system if the user intends to resample an image with gregister or geotran, or
transform coordinates from the reference coordinate system to the input
image coordinate system. 

.nf
    xin = f (xref, yref)
    yin = g (xref, yref)
.fi

If on the other hand the user wishes to transform coordinates from the
input image coordinate system to the reference coordinate system then he or she
must reverse the roles of the reference and input coordinates as defined above,
and compute the inverse transformation.


The functions f and g are either a power series polynomial or a Legendre or
Chebyshev polynomial surface of order \fIxxorder\fR and \fIxyorder\fR in x
and \fIyxorder\fR and \fIyyorder\fR in y.

Several polynomial cross terms options are available. Options "none",
"half", and "full" are illustrated below for a quadratic polynomial in
x and y.

.nf
xxterms = "none", xyterms = "none"
xxorder = 3, xyorder = 3, yxorder = 3, yyorder = 3

   xin = a11 + a21 * xref + a12 * yref +
         a31 * xref ** 2 + a13 * yref ** 2
   yin = a11' + a21' * xref + a12' * yref +
         a31' * xref ** 2 + a13' * yref ** 2

xxterms = "half", xyterms = "half"
xxorder = 3, xyorder = 3, yxorder = 3, yyorder = 3

   xin = a11 + a21 * xref + a12 * yref +
         a31 * xref ** 2 + a22 * xref * yref + a13 * yref ** 2
   yin = a11' + a21' * xref + a12' * yref +
         a31' * xref ** 2 + a22' * xref * yref + a13' * yref ** 2

xxterms = "full", xyterms = "full"
xxorder = 3, xyorder = 3, yxorder = 3, yyorder = 3

   xin = a11 + a21 * xref + a31 * xref ** 2 +
         a12 * yref + a22 * xref * yref +  a32 * xref ** 2 * yref +
         a13 * yref ** 2 + a23 * xref *  yref ** 2 +
         a33 * xref ** 2 * yref ** 2
   yin = a11' + a21' * xref + a31' * xref ** 2 +
         a12' * yref + a22' * xref * yref +  a32' * xref ** 2 * yref +
         a13' * yref ** 2 + a23' * xref *  yref ** 2 +
         a33' * xref ** 2 * yref ** 2
.fi

If the \fBfitgeometry\fR parameter is anything other than "general", the  order
parameters assume the value 2 and the cross terms switches assume the value
"none", regardless of the values set by the user. The computation can be done in
either real or double precision by setting \fIcalctype\fR. Automatic pixel
rejection may be enabled by setting \fmaxiter\fR > 0 and \fIreject\fR to some
number greater than 0.

\fIXmin\fR, \fIxmax\fR, \fIymin\fR and \fIymax\fR define the region of
validity of the fit in the reference coordinate system and must be set by
the user. These parameters can be used to reject out of range data before the
actual fitting is done.

GEOMAP may be run interactively by setting \fIinteractive\fR = yes and
inputting commands by the use of simple keystrokes.
In interactive mode the user has the option of changing the
fit parameters and displaying the data graphically until a satisfactory
fit has been achieved. The available keystroke commands are listed
below.

.nf
?	Print options
f	Fit the data and graph with the current graph type (g, x, r, y, s)
g	Graph the data and the current fit
x,r	Graph the x fit residuals versus x and y respectively
y,s	Graph the y fit residuals versus x and y respectively
d,u	Delete or undelete the data point nearest the cursor
o	Overplot the next graph
c	Toggle the constant x, y plotting option
t       Plot a line of constant x, y through the nearest data point	
l	Print xshift, yshift, xmag, ymag, xrotate, yrotate
q	Exit the interactive curve fitting
.fi

The parameters listed below can be changed interactively with simple colon
commands. Typing the parameter name alone will list the current value.

.nf
:show				List parameters
:fitgeometry			Fitting geometry (shift,xyscale,rotate,
				rscale,rxyscale,general)
:function [value]	        Fitting function (chebyshev,legendre,
                                polynomial)
:xxorder :xyorder [value]	X fitting function xorder, yorder
:yxorder :yyorder [value]	Y fitting function xorder, yorder
:xxterms :yxterms [n/h/f]	X, Y fit cross terms type
:maxiter [value]		Maximum number of rejection iterations
:reject [value]			Rejection threshold
.fi

The final fit is stored in a simple text file in a format suitable for use
by the GREGISTER or GEOTRAN tasks.

If \fIverbose\fR  is "yes", various pieces of useful information are printed
to the terminal as the task proceeds. If \fIresults\fR is set to a file name
then the input coordinates, the fitted coordinates, and the residuals of
the fit are written to that file.

The transformation computed by the "general" fitting geometry is arbitrary
and does not correspond to a physically meaningful model. However the computed
coefficients for the linear term can be given a simple geometrical geometric
interpretation for all the fitting geometries as shown below.

.nf
	fitting geometry = general (linear term)
	    xin = a + b * xref + c * yref
	    yin = d + e * xref + f * yref

	fitting geometry = shift
	    xin = a + xref
	    yin = d + yref

	fitting geometry = xyscale
	    xin = a + b * xref
	    yin = d + f * yref

	fitting geometry = rotate
	    xin = a + b * xref + c * yref
	    yin = d + e * xref + f * yref
	    b * f - c * e = +/-1
	    b = f, c = -e or b = -f, c = e

	fitting geometry = rscale
	    xin = a + b * xref + c * yref
	    yin = d + e * xref + f * yref
	    b * f - c * e = +/- const
	    b = f, c = -e or b = -f, c = e

	fitting geometry = rxyscale
	    xin = a + b * xref + c * yref
	    yin = d + e * xref + f * yref
	    b * f - c * e = +/- const
.fi

The coefficients can be interpreted as follows. Xref0, yref0, xin0, yin0
are the origins in the reference and input frames respectively. Orientation
and skew are the rotation of the x and y axes and their deviation from
perpendicularity respectively. Xmag and ymag are the scaling factors in x and
y and are assumed to be positive.

.nf
	general (linear term)
	    xrotation = rotation - skew / 2
	    yrotation = rotation + skew / 2
	    b = xmag * cos (xrotation)
	    c = ymag * sin (yrotation)
	    e = -xmag * sin (xrotation)
	    f = ymag * cos (yrotation)
	    a = xin0 - b * xref0 - c * yref0 = xshift
	    d = yin0 - e * xref0 - f * yref0 = yshift

	shift
	    xrotation = 0.0,  yrotation = 0.0
	    xmag = ymag = 1.0
	    b = 1.0
	    c = 0.0
	    e = 0.0
	    f = 1.0
	    a = xin0 - xref0 = xshift
	    d = yin0 - yref0 = yshift

	xyscale
	    xrotation 0.0 / 180.0 yrotation = 0.0
	    b = + /- xmag
	    c = 0.0
	    e = 0.0
	    f = ymag
	    a = xin0 - b * xref0 = xshift
	    d = yin0 - f * yref0 = yshift

	rscale
	    xrotation = rotation + 0 / 180, yrotation = rotation
	    mag = xmag = ymag
	    const = mag * mag
	    b = mag * cos (xrotation)
	    c = mag * sin (yrotation)
	    e = -mag * sin (xrotation)
	    f = mag * cos (yrotation)
	    a = xin0 - b * xref0 - c * yref0 = xshift
	    d = yin0 - e * xref0 - f * yref0 = yshift

	rxyscale
	    xrotation = rotation + 0 / 180, yrotation = rotation
	    const = xmag * ymag
	    b = xmag * cos (xrotation)
	    c = ymag * sin (yrotation)
	    e = -xmag * sin (xrotation)
	    f = ymag * cos (yrotation)
	    a = xin0 - b * xref0 - c * yref0 = xshift
	    d = yin0 - e * xref0 - f * yref0 = yshift
.fi


.ih
EXAMPLES
1. Compute the linear transformation between coordinate systems.
   A record called "m51.coo" will be written in the database
   file "database".


.nf
	cl> geomap m51.coo database 1. 512. 1. 512.
.fi

2. Compute the 3rd order transformation in x and y between two
   coordinate systems.  A record called "m51.coo" will be written in
   the database file "database". This record supersedes the one
   of the same name written in example 1.

.nf
	cl> geomap m51.coo database 1. 512. 1. 512. xxo=4 xyo=4 \
	>>> yxo=4 yyo=4 xxt=full yxt=full inter-
.fi

3. Register a 500 by 500 image of m51 to an 800 by 800 image of the same
field taken with a different instrument, and display the original
800 by 800 image and the transformed image. Use the default fitting parameters.

.nf
	cl> geomap m51.coo database 1.0 800.0 1.0 800.0
	cl> gregister m51.500 m51.500.out database m51.coo
	cl> display m51.800 1 fi+
	cl> display m51.500.out 2 fi+
.fi

4. Use the above transform to transform a list of object pixel coordinates
in the m51.800 image to their pixel coordinates in the m51.500 system.

.nf
	cl> geoxytran m51.800.xy m51.500.xy database m51.coo
.fi

5. Transform object pixel coordinates in the m51.500 image to their
pixel coordinates in the m51.800 image. Note that to do this the roles
of the reference and input coordinates defined in example 3 must be
reversed and the inverse transform must be computed.

.nf
	cl> fields m51.coo 3,4,1,2 > m51.coo.inv
	cl> geomap m51.coo.inv database 1.0 512.0 1.0 512.0
	cl> geoxytran m51.512.xy m51.800.xy database m51.coo.inv
.fi

6. Compute 3 different transforms, store them in the same database file,
and use them to transform 3 different images.  Use the original image names as
the database record names.

.nf
	cl> geomap coo1,coo2,coo3 database 1. 512. 1. 512. \
	>>> transforms=im1,im2,im3
	cl> gregister im1,im2,im3  im1.out,im2.out,im3.out database \
	>>> im1,im2,im3
.fi

.ih
BUGS

The user should be aware that for high order fits the "polynomial" basis
functions become very unstable. Switching to the "legendre" or "chebyshev"
polynomials and/or going to double precision will usually cure the problem.

.ih
SEE ALSO
imshift, magnify, rotate, imlintran, gregister, geotran, geoxytran
.endhelp