| wcsmap | images.immatch | wcsmap |
wcsmap -- compute the spatial transformation function required to register a list of images using WCS information
wcsmap input reference database
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 set to "none", regardless of the values of the xxorder , xyorder , xxterms , yxorder , yyorder and yxterms parameters set by the user.
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 is set to "half" or "full" and fitgeometry is "general" then a distortion surface is fit to the residuals from the linear portion of the fit.
WCSMAP computes the spatial transformation function required to map the coordinate system of the reference image reference to the coordinate system of the input image input , and stores the computed function in the output text database file database . The input and reference images must be one- or two-dimensional and must have the same dimensionality. The input image and output text database file can be input to the REGISTER or GEOTRAN tasks to perform the actual image registration. WCSMAP assumes that the world coordinate systems in the input and reference image headers are accurate and that the two systems are compatible, e.g. both images have the same epoch sky projection world coordinate systems or both are spectra whose coordinates are in the same units.
WCSMAP computes the required spatial transformation by matching the logical x and y pixel coordinates of a grid of points in the input image with the logical x and y pixels coordinates of the same grid of points in the reference image, using world coordinate information stored in the two image headers. The coordinate grid consists of nx * ny points evenly distributed over the logical pixel space of interest in the reference image defined by the xmin , xmax , ymin , ymax parameters. The logical x and y pixel reference image coordinates are transformed to the reference image world coordinate system defined by wcs , using the wcs information in the reference image header. The reference image world coordinates are then transformed to logical x and y pixel coordinates in the input image, using world coordinate system information stored in the input image header.
The computed reference and input logical coordinates and the world coordinates are written to a temporary output coordinates file which is deleted on task termination. The x and y coordinates are written using the xformat and yformat and the wxformat and wxformat parameters respectively. If these formats are undefined and, in the case of the world coordinates a format attribute cannot be read from either the reference or the input images, the coordinates are output in %g format with min_sigdigits digits of precision. If the reference and input images are 1D then all the output logical and world y coordinates are set to 1.
WCSMAP computes a spatial transformation of the following form.
xin = f (xref, yref)
yin = g (xref, yref)
The functions f and g are either a power series polynomial or a Legendre or Chebyshev polynomial surface of order xxorder and xyorder in x and yxorder and yyorder in y. Cross terms are optional.
Several polynomial cross terms options are avaible. Options "none", "half", and "full" are illustrated below for a quadratic polynomial in x and y.
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
If the fitgeometry 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 the calctype parameter. Automatic pixel rejection may be enabled by setting the reject parameter to some number > 0.0.
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.
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
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 orientation 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.
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
Xmin , xmax , ymin and ymax define the region of validity of the fit as well as the limits of the grid in the reference coordinate system and must be set by the user. These parameters are used to reject out of range data before the actual fitting is done.
Each computed transformation is written to the output file database in a record whose name is either specified by the user via the transforms parameter or defaults the name of the input image. The database file is opened in append mode and new records are written to the end of the existing file. If more that one record of the same name is written to the database file, the last record written is the valid record, i.e. the one that will be used by the REGISTER or GEOTRAN tasks.
WCSMAP will terminate with an error if the reference and input images are not both either 1D or 2D. If the world coordinate system information cannot be read from either the reference or input image header, the requested transformations from the world <-> logical coordinate systems cannot be compiled for either or both images, or the world coordinate systems of the reference and input images are fundamentally incompatible in some way, the output logical reference and input image coordinates are both set to a grid of points spanning the logical pixel space of the input, not the reference image. This grid of points defines an identity transformation which will leave the input image unchanged if applied by the REGISTER or GEOTRAN tasks.
If verbose is "yes" then messages about the progress of the task as well as warning messages indicating potential problems are written to the standard output. If results is set to a file name then the input coordinates, the fitted coordinates, and the residuals of the fit are written to that file.
WCSMAP may be run interactively by setting the interactive parameter to "yes". In interactive mode the user has the option of viewing the fit, changing the fit parameters, deleting and undeleting points, and replotting the data until a satisfactory fit has been achieved.
In interactive mode the the following cursor commands are currently available.
Interactive Keystroke Commands
? 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
The parameters listed below can be changed interactively with simple colon commands. Typing the parameter name alone will list the current value.
Colon Parameter Editing Commands
: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
:reject [value] Rejection threshold
A format specification has the form "%w.dCn", where w is the field width, d is the number of decimal places or the number of digits of precision, C is the format code, and n is radix character for format code "r" only. The w and d fields are optional. The format codes C are as follows:
b boolean (YES or NO)
c single character (c or '\c' or '\0nnn')
d decimal integer
e exponential format (D specifies the precision)
f fixed format (D specifies the number of decimal places)
g general format (D specifies the precision)
h hms format (hh:mm:ss.ss, D = no. decimal places)
m minutes, seconds (or hours, minutes) (mm:ss.ss)
o octal integer
rN convert integer in any radix N
s string (D field specifies max chars to print)
t advance To column given as field W
u unsigned decimal integer
w output the number of spaces given by field W
x hexadecimal integer
z complex format (r,r) (D = precision)
Conventions for w (field width) specification:
W = n right justify in field of N characters, blank fill
-n left justify in field of N characters, blank fill
0n zero fill at left (only if right justified)
absent, 0 use as much space as needed (D field sets precision)
Escape sequences (e.g. "\n" for newline):
\b backspace (not implemented)
\f formfeed
\n newline (crlf)
\r carriage return
\t tab
\" string delimiter character
\' character constant delimiter character
\\ backslash character
\nnn octal value of character
Examples
%s format a string using as much space as required
%-10s left justify a string in a field of 10 characters
%-10.10s left justify and truncate a string in a field of 10 characters
%10s right justify a string in a field of 10 characters
%10.10s right justify and truncate a string in a field of 10 characters
%7.3f print a real number right justified in floating point format
%-7.3f same as above but left justified
%15.7e print a real number right justified in exponential format
%-15.7e same as above but left justified
%12.5g print a real number right justified in general format
%-12.5g same as above but left justified
%h format as nn:nn:nn.n
%15h right justify nn:nn:nn.n in field of 15 characters
%-15h left justify nn:nn:nn.n in a field of 15 characters
%12.2h right justify nn:nn:nn.nn
%-12.2h left justify nn:nn:nn.nn
%H / by 15 and format as nn:nn:nn.n
%15H / by 15 and right justify nn:nn:nn.n in field of 15 characters
%-15H / by 15 and left justify nn:nn:nn.n in field of 15 characters
%12.2H / by 15 and right justify nn:nn:nn.nn
%-12.2H / by 15 and left justify nn:nn:nn.nn
\n insert a newline
Additional information on IRAF world coordinate systems including more detailed descriptions of the "logical", "physical", and "world" coordinate systems can be found in the help pages for the WCSEDIT and WCRESET tasks. Detailed documentation for the IRAF world coordinate system interface MWCS can be found in the file "iraf$sys/mwcs/MWCS.hlp". This file can be formatted and printed with the command "help iraf$sys/mwcs/MWCS.hlp fi+ | lprint". Information on the spectral coordinates systems and their suitability for use with WCSXYMATCH can be obtained by typing "help specwcs | lprint". Details of the FITS header world coordinate system interface can be found in the document "World Coordinate Systems Representations Within the FITS Format" by Hanisch and Wells, available from our anonymous ftp archive.
1. Compute the spatial transformation required to match a radio image to an X-ray image of the same field field using a 100 point coordinate grid and a simple linear transformation. Both images have accurate sky projection world coordinate systems. Print the output world coordinates in the coords file in hh:mm:ss.ss and dd:mm:ss.s format. Run geotran on the results to do the actual registration.
cl> wcsmap radio xray geodb wxformat=%12.2H wyformat=%12.1h \ interactive- cl> geotran radio radio.tran geodb radio
2. Repeat the previous command but begin with a higher order fit and run the task in interactive mode in order to examine the fit residuals.
cl> wcsmap radio xray geodb wxformat=%12.2H wyformat=%12.1h \
xxo=4 xyo=4 xxt=half yxo=4 yyo=4 yxt=half
... a plot of the fit appears
... type x and r to examine the residuals of the x
surface fit versus x and y
... type y and s to examine the residuals of the y
surface fit versus x and y
... delete 2 deviant points with the d key and
recompute the fit with the f key
... type q to quit and save the fit
cl> geotran radio radio.tran geodb radio
3. Repeat example 1 but assign a user name to the transform.
cl> wcsmap radio xray geodb transforms="m82" wxformat=%12.2H \ wyformat=%12.1h interactive- cl> geotran radio radio.tran geodb m82
wcstran,xregister,wcsxymatch,geomap,register,geotran