.help psfmatch Oct94 images.immatch
.ih
NAME
psfmatch -- match the point spread functions of 1 and 2D images
.ih
USAGE
psfmatch input reference psfdata kernel 
.ih
PARAMETERS
.ls input
The list of input images to be matched.
.le
.ls reference
The list of reference images to which the input images are to be matched if
\fIconvolution\fR = "image", or the list of reference image psfs if 
\fIconvolution\fR = "psf". The reference image psf must be broader than the
input image psf in at least one dimension.
The number of reference images/psfs must be one or equal to the number of
input images.
.le
.ls psfdata
The list of objects used to compute the psf matching function if
\fIconvolution\fR is "image", or the list of input image psfs if 
\fIconvolution\fR is "psf". In the former case \fIpsfdata\fR may be:
1) a string containing the x and y coordinates of a single object,
e.g. "51.0 105.0" or 2) the name of a text file containing a list of
objects, and the number of objects
files must equal the number of reference images. In the latter case
the number of input psf images must equal the number of input images.
.le
.ls kernel
The list of input/output psf matching function images to be convolved with the
input images to produce the output images. The number of kernel images
must equal the number of input images.
.le
.ls output = ""
The list of output matched images. If \fIoutput\fR is the NULL string
then the psf matching function is computed for each input image and written to
\fIkernel\fR but no output images are written. If \fIoutput\fR is not NULL
then the number of output images must equal the number of input images.
.le
.ls convolution = "image"
The algorithm used to compute the psf matching function. The options are:
.ls image
The psf matching function is computed directly from the reference and input
image data using the objects specified in \fIpsfdata\fR, the data
regions specified by \fIdnx\fR, \fIdny\fR, \fIpnx\fR, and \fIpny\fR,
and the convolution theorem.
.le
.ls psf   
The psf matching function is computed directly from pre-computed
reference and input image psfs using the convolution theorem.
.le
.ls kernel
No psf matching function is computed. Instead the psf matching function
is  read from the input image \fIkernel\fR.
.le
.le
.ls dnx = 31, ls dny = 31
The x and y width of the data region to be extracted around each object. The
data region should be big enough to include both object and sky data.
\fIDnx\fR and \fIdny\fR are not used if \fIconvolution\fR is "psf" or
"kernel".
.le
.ls pnx = 15, pny = 15
The x and y width of the psf matching function to be computed which must be
less than \fIdnx\fR and \fIdny\fR respectively. The psf
matching function should be kept as small as possible to minimize
the time required to compute the output image.
\fIPnx\fR and \fIPny\fR are not used if \fIconvolution\fR is "psf" or
"kernel".
.le
.ls center = yes
Center the objects in \fIpsfdata\fR before extracting the data from the
input and reference images. Centering should be turned off if the objects
are non-stellar and do not have well-defined centers.
Centering is turned off if \fIconvolution\fR is "psf" or
"kernel".
.le
.ls background = median
The default background function to be subtracted from the input
and reference image data in each object region before the
psf matching function is computed. The background is computed using
data inside the data extraction region defined by \fIdnx\fR and \fIdny\fR
but outside the kernel region defined by \fIpnx\fR and \fIpny\fR.
Background fitting is turned off if \fIconvolution\fR is "psf" or
"kernel".
The options are:
.ls none
no background subtraction is done.
.le
.ls "insky refsky"
the numerical values of insky and refsky are subtracted from the
input and reference image respectively.
.le
.ls mean
the mean of the input and reference image region is computed and subtracted
from the image data.
.le
.ls median
the median of the input and reference image region is computed and subtracted
from the data.
.le
.ls plane
a plane is fit to the input and reference image region and subtracted
from the data.
.le
.le
.ls loreject = INDEF, ls hireject = INDEF
The k-sigma rejection limits for removing the effects of bad data from the
background fit.
.le
.ls apodize = 0.0
The fraction of the input and reference image data endpoints in x and y
to apodize with a
cosine bell function before the psf matching function is computed.
Apodizing is turned off if \fIconvolution\fR is "psf" or
"kernel".
.le
.ls fluxratio = INDEF
The ratio of the integrated flux of the reference objects to the integrated
flux of the input objects.
By default \fIfluxratio\fR is computed directly from the input data.
.le
.ls filter = "replace"
The filter used to remove high frequency noise from the psf
matching function. Filtering is not performed if \fIconvolution\fR
is "kernel". The options are:
.ls cosbell
apply a cosine bell taper to the psf matching function in frequency space. 
.le
.ls replace
replace the high-frequency low signal-to-noise components of the psf matching
function with a gaussian model computed from the low frequency
high signal-to-noise components of the matching function.
.le
.ls model
replace the entire psf matching function with a gaussian model fit to the
low frequency high signal-to-noise components of the matching function.
.le
.le
.ls sx1 = INDEF, sx2 = INDEF, sy1 = INDEF, sy2 = INDEF
The limits of the cosine bell taper in frequency space. Frequency components
inside sx1 and sy1 are unaltered. Frequency components outside sx2 and sy2
are set to 0.0. By default sx1 and sy1 are set to 0.0,
and sx2 and sy2 are set to the largest frequency present in the data.
.le
.ls radsym = no
Compute a radially symmetric cosine bell function ?
.le
.ls threshold = 0.2
The low frequency cutoff in fraction of the total input image spectrum
power for the filtering options "replace" and "model".
.le
.ls normfactor = 1.0
The total power in the computed psf matching function \fIkernel\fR. By default
the psf matching function is normalized.  If \fInormfactor\fR
is set to INDEF, then the total power is set to \fIfluxratio\fR.
\fINormfactor\fR is not used if \fIconvolution\fR is set "kernel".
.le
.ls boundary_type = "nearest"
The boundary extension algorithm used to compute the output matched
image.  The options are:
.ls nearest
use the value of the nearest boundary pixel.
.le
.ls constant
use a constant value.
.le
.ls reflect
generate a value by reflecting about the boundary.
.le
.ls wrap
generate a value by wrapping around to the opposite side of the image.
.le
.le
.ls constant = 0.0
The default constant for constant boundary extension.
.le
.ls interactive = no
Compute the psf matching function for each image
interactively using graphics cursor and, optionally, image cursor input.
.le
.ls verbose
Print messages about the progress of the task in non-interactive mode.
.le
.ls graphics = "stdgraph"
The default graphics device.
.le
.ls display = "stdimage"
The default image display device.
.le
.ls gcommands = ""
The default graphics cursor.
.le
.ls icommands = ""
The default image display cursor.
.le

.ih
DESCRIPTION

PSFMATCH computes the convolution kernel required to match the
point-spread functions
of the input images \fIinput\fR to the point-spread functions of
the reference images \fIreference\fR using either the image data 
or pre-computed psfs and the convolution theorem.
The computed psf matching functions are stored in the \fIkernel\fR images.
If a non-NULL list of output images \fIoutput\fR is
specified the input images are
convolved with the kernel images to produce a list of psf matched output
images. PSFMATCH requires
that the input and reference images be spatially registered
and that the reference images have poorer resolution (broader PSF)
than the input images in at least one dimension.

If \fIconvolution\fR = "image", the matching function is computed directly
from the input and reference image data using the objects listed in
\fIpsfdata\fR and the convolution theorem as described in the ALGORITHMS
section. \fIpsfdata\fR is interpreted as either: 1) a
string defining the coordinates of a single object e.g. "103.3 189.2" or 2)
the name of a text file containing the coordinates of one or 
more objects, one object per line, with the x and y coordinates
in columns 1 and 2 respectively.  The object coordinates, the
size of the data region to be extracted \fIdnx\fR
by \fIdny\fR, and the size of the kernel to be computed \fIpnx\fR and
\fIpny\fR, determine 
the input and reference image regions used to compute the psf matching
function.
These image regions should be selected with care. Ideal regions 
contain a single high signal-to-noise unsaturated star which has no close
neighbors and is well centered on a pixel.

If \fIcenter\fR is "yes" and \fIconvolution\fR is "image", the objects
in \fIpsfdata\fR are centered before
the data region is extracted.  Centering should be on if the objects
are stellar, particularly if their coordinates were read from the image
display cursor. Centering should be off if the objects are non-stellar and
do not have well-defined centers.

If the \fIbackground\fR fitting algorithm is other than "none" and
\fIconvolution\fR is "image", the background for each object is fit using 
data inside the region defined by
\fIdnx\fR and \fIdny\fR but outside the region defined by
\fIpnx\fR by \fIpny\fR. Bad data can be removed from the
background fit by setting the parameters \fIloreject\fR and \fIhireject\fR.
A cosine bell function is applied to the edges of the data region
after background fitting but before computing the psf matching function
if the \fIapodize\fR parameter is > 0.0.

If \fIpsfdata\fR contains more than one object, the extracted image data
is weighted by the total intensity in the extracted region after
background subtraction, and averaged to produce a single smoothed
data region for each reference and input image.

If \fIconvolution\fR = "psf",
the psf matching function is computed directly from the input image
and reference
image point-spread functions
using the convolution theorem as described in the ALGORITHMS section.
In this case  \fIpsfdata\fR is the list of input image psfs  and
\fIreference\fR are the corresponding reference image psfs written by
by some external psf modeling task. 
If \fIconvolution\fR is "psf",
centering and background fitting
are assumed to have been performed by the psf modeling task and are not
performed by PSFMATCH.

PSFMATCH requires that the total power in the psf matching function
before normalization be the ratio
of the integrated flux of the reference image/psf over the integrated
flux of the input image/psf. If \fIfluxratio\fR is INDEF, PSFMATCH
estimates this number internally as described in the ALGORITHMS section,
otherwise the \fIfluxratio\fR is set to the value supplied by the user.

If \fIconvolution\fR is "kernel", PSFMATCH reads the psf matching function
from the images in \fIkernel\fR  which were either
created during a previous run of PSFMATCH or by a separate task.

PSFMATCH provides several options for filtering out the ill-behaved
noise-dominated high frequency components of the psf matching function
that are produced when the ratio of reference / input image of psf
fourier transforms is taken.

If \fIfilter\fR is set to "cosbell", a cosine bell function
with a taper defined by \fIsx1\fR, \fIsx2\fR, \fIsy1\fR, and \fIsy2\fR and
symmetry defined by \fRradsym\fR is applied to
the psf matching function in frequency space. This filter
sets all the frequency components greater than \fIsx2\fR and \fIsy2\fR
to 0.0 and leaves all frequency components inside \fIsx1\fR and \fIsy1\fR
unaltered. Users should exercise this option with caution as the effect
of the filtering process can be to significantly
broaden the computed psf matching function as described in the ALGORITHMS
section.

An alternative approach to dealing with the noisy
high frequency components of the psf
matching function it is to replace them with a reasonable guess. If the
matching function is approximately gaussian then its fourier transform is also
approximately gaussian and the low frequency components can be modeled
reliably with an elliptical gaussian function. The model derived from the low
frequency components of the matching can then be used to replace the high
frequency components.
If \fIfilter\fR is set to "replace", those high frequency components
of the matching function  which have less than a fraction
\fIthreshold\fR of their total power in the equivalent high frequency
components of the divisor or input image transform,
are replaced by a model computed by fitting a gaussian to the low frequency
components of the matching function, as described in the ALGORITHMS section.
If \fIfilter\fR = "model" then the entire psf matching function
is replaced with the best fitting gaussian model.

Another problem can arise during the computation of the psf matching
function . Occasionally it is not possible by means of a single execution
of PSFMATCH to match the reference and input image psfs. An example
of this situation
is the case where the seeing of the reference and input images
was comparable but the declination guiding error in the reference
image was larger than the error in the input image.
In this case input image  needs to be convolved to the resolution of 
the reference image. However it is also the case
that the guiding error in ra in the input image is greater than the guiding
error  in ra in the reference image. In this case the reference image needs
to be convolved to the resolution of the input image along the other axis.
If no corrective action is taken by the task, the 
first time PSFMATCH is run the values of the psf matching function along
the ra axis will be greater than the computed fluxratio, resulting in
unrealistic action
along this axis. PSFMATCH avoids this situation by internally limiting
the psf matching function to a maximum value of fluxratio computed as described
above. 

By default the psf matching function is normalized to unit power before 
output. This may not be what is desired since if carefully computed the
internally computed quantity a contains information about differences
in exposure time, transparency, etc. If \fInormfactor\fR is set to
a number of INDEF, the total power of the psf matching function will be
set to that value of \fIfluxratio\fR respectively.

If a list of output images names has been supplied then the computed
psf matching function is applied to the input images to produce
the output images using the boundary extension algorithm
defined by \fIboundary\fR and \fIconstant\fR.

In non-interactive mode the parameters are set at task startup time and
the input images are processed sequentially. If the \fIverbose\fR flag
is set messages about the progress of the task are printed on he 
screen as the task is running.

In interactive mode the user can mark the regions to be used to compute
the psf matching function on the image display, show/set the data
and algorithm parameters, compute, recompute, and plot the psf matching
function and its accompanying fourier spectrum, and experiment with the
various filtering and modeling options.

.ih
CURSOR COMMANDS

The following graphics cursor commands are currently available in
PSFMATCH.

.nf
	Interactive Keystroke Commands


?	Print help 
:	Colon commands
k	Draw a contour plot of the psf matching kernel
p	Draw a contour plot of the psf matching kernel fourier spectrum
x	Draw a column plot of the psf matching kernel / fourier spectrum
y	Draw a line plot of the psf matching kernel / fourier spectrum
r	Redraw the current plot
f	Recompute the psf matching kernel
w	Update the task parameters
q	Exit


	Colon Commands


:mark	[file]		Mark objects on the display
:show			Show current values of the parameters


	Show/Set Parameters


:input	    [string]	    Show/set the current input image name
:reference  [string]	    Show/set the current reference image/psf name
:psf	    [file/string]   Show/set the objects/input psf list
:psfimage   [string]	    Show/set the current input psf name
:kernel	    [string]	    Show/set the current psf matching kernel name
:output     [string]	    Show/set the current output image name

:dnx	    [value]	    Show/set x width of data region(s) to extract
:dny	    [value]	    Show/set y width of data region(s) to extract
:pnx	    [value]	    Show/set x width of psf matching kernel
:pny	    [value]	    Show/set y width of psf matching kernel
:center	    [yes/no]	    Show/set the centering switch
:background [string]        Show/set the background fitting function
:loreject   [value]	    Show/set low side k-sigma rejection parameter
:hireject   [value]	    Show/set high side k-sigma rejection parameter
:apodize    [value]	    Show/set percent of endpoints to apodize

:filter	    [string]	    Show/set the filtering algorithm
:fluxratio  [value]	    Show/set the reference/input psf flux ratio
:sx1	    [value]	    Show/set inner x frequency for cosbell filter
:sx2	    [value]	    Show/set outer x frequency for cosbell filter
:sy1	    [value]	    Show/set inner y frequency for cosbell filter
:sy2	    [value]	    Show/set outer y frequency for cosbell filter
:radsym	    [yes/no]        Show/set radial symmetry for cosbell filter
:threshold  [value]	    Show/set %threshold for replace/modeling filter
:normfactor [value]	    Show/set the kernel normalization factor
.fi

.ih
ALGORITHMS

The problem of computing the psf matching function can expressed
via the convolution theorem as shown below.
In the following expressions r is the reference
image data or reference image psf, i is the input image data or input image
psf, k is the unit power psf matching
function,
a is a scale factor specifying the ratio of the total
power in the reference data or psf to the total power in the input data or
psf, * is the convolution operator, and FT is the fourier transform operator.

.nf
	r = ak * d
	R = FT (r)
	I = FT (i)
	aK = R / I
	ak = FT (aK)
.fi

The quantity ak is the desired psf matching function and aK is its fourier
transform.

If the background was accurately removed from the image or psf data before the
psf matching function was computed, the quantity a is simply the central
frequency component of the computed psf matching function aK as shown below.

.nf
	aK[0,0] = a = sum(r) / sum(i)
.fi

If the background was not removed from the image or psf data before the
psf matching function was computed the previous expression is not valid.
The computed aK[0,0] will include an offset and a must be estimated
in some other manner. The approach taken by PSFMATCH in this circumstance
is to fit a gaussian model to the absolute value of 1st and 2nd frequencies
of R and I along the x and y axes independently, average the fitted x and y
amplitudes, and set aK[0,0] to the ratio of the resulting fitted amplitudes
as shown below.

.nf
	      a = amplitude (R) / amplitude (I)
	        = (sum(r) - sum(skyr)) / (sum(i) - sum(skyi))  
	      aK[0,0] = a
.fi

This approach will work well as long as the image data or psf is reasonably
gaussian but may not work well in arbitrary image regions. If the user is
dissatisfied with either of the techniques described above they can
set aK[0,0] to a pre-determined value of their own.

If a filter is applied to the computed psf matching function in frequency
space then instead of computing

.nf
	       ak = FT (aK)
.fi

PSFMATCH actually computes

.nf
	       ak' = FT (aKF) = ak * f
.fi

where F is the applied filter in frequency space and f is its
fourier transform. Care should be taken in applying any filter.
For example if F is the step function, then ak' will be the desired kernel
ak convolved with f, a sinc function of frequency 2 * PI / hwidth where
hwidth is the half-width of the step function, and the resulting k'
will be too broad.

If the user chooses to replace the high frequency components of the psf
matching function with a best guess, PSFMATCH performs the following
steps:

.nf
1) fits an elliptical gaussian to those frequency components of the fourier
spectrum of aK for which for which the amplitude of I is greater
than threshold * I[0,0] to determine the geometry of the ellipse

2) uses the fourier shift theorem to preserve the phase information in the
model and solve for any x and y shifts

3) replace those frequency components of aK for which the fourier spectrum
of I is less than threshold * I[0,0] with the model values

		or alternatively

replace all of aK with the model values
.fi

.ih
EXAMPLES

1. Psf match a list of input images taken at different epochs with variable
seeing conditions to a reference image with the poorest seeing by marking
several high signal-to-noise isolated stars on the displayed reference image
and computing the psf matching function directly from the input and reference
image data. User makes two runs with psfmatch one to compute and check the
kernel images and one to match the images.

.nf
	cl> display refimage 1 fi+

	cl> rimcursor > objects

	cl> psfmatch @inimlist refimage objects @kernels dnx=31 \
	    dny=31 pnx=15 pny=15

	cl> imstat @kernels

	cl> psfmatch @inlist refimage objects @kernels          \
	    output=@outlist convolution="kernel"
.fi

2. Psf match two spectra using a high signal-to-noise portion of the
data in the middle of the spectrum. Since the spectra are registered
spatially and there is little data available for background fitting the
user chooses to turn centering off and set the backgrounds manually.

.nf
	cl> psfmatch inspec refspec "303.0 1.0" kernel         \
	    output=outspec dnx=31 dny=31 pnx=15 pny=15 center- \
	    back="403.6 452.0"
.fi

3. Psf match two images using psf functions inpsf and refpsf computed with
the daophot package phot/psf/seepsf tasks. Since the kernel is fairly
large use the stsdas fourier package task fconvolve to do the actual
convolution. The boundary extension algorithm in fconvolve is equivalent
to setting the psfmatch boundary extension parameters boundary and
constant to "constant" and "0.0" respectively.

.nf
	cl> psfmatch inimage refpsf inpsf kernel convolution=psf

	cl> fconvolve inimage kernel outimage
.fi

4. Psf match two images interactively using the image data itself to
compute the psf matching function.

.nf
	cl> psfmatch inimage refimage objects kernel interactive+

	    ... a contour plot of the psf matching function appears
		with the graphics cursor ready to accept commands

            ... type x and y to get line and column plots of the psf
                matching function at various points and k to return
                to the default contour plot

	    ... type ? to get a list of the available commands

	    ... type :mark to define a new set of objects

	    ... type f to recompute the psf matching function using
                the new objects

 	    ... increase the data window to 63 pixels in x and y
                with the :dnx 63 and :dny 63 commands, at the
                same time increase the psf function size to 31 with
		the colon commands :pnx 31 and :pny 31

	    ... type f to recompute the psf matching function using
                the new data and kernel windows

	    ... type q to quit the task, and q again to verify the previous
                q command
.fi

.ih
TIME REQUIREMENTS
.ih
BUGS
.ih
SEE ALSO
convolve, gauss, stsdas.fconvolve, digiphot.daophot.psf
.endhelp
