This is a set of functions for resampling gridded data (e.g. an image) under the control of a geometrical transformation, which is specified by a Mapping. The functions operate on a pair of data grids (input and output), each of which may have any number of dimensions. Resampling may be restricted to a specified region of the output grid. An associated grid of error estimates associated with the input data may also be supplied (in the form of variance values), so as to produce error estimates for the resampled output data. Propagation of missing data (bad pixels) is supported.

You should use a resampling function which matches the numerical type of the data you are processing by replacing $<$X$>$ in the generic function name AST_RESAMPLE$<$X$>$ by an appropriate 1- or 2-character type code. For example, if you are resampling data with type REAL, you should use the function AST_RESAMPLER (see the " Data Type Codes" section below for the codes appropriate to other numerical types).

Resampling of the grid of input data is performed by transforming the coordinates of the centre of each output grid element (or pixel) into the coordinate system of the input grid. Since the resulting coordinates will not, in general, coincide with the centre of an input pixel, sub-pixel interpolation is performed between the neighbouring input pixels. This produces a resampled value which is then assigned to the output pixel. A choice of sub-pixel interpolation schemes is provided, but you may also implement your own.

This algorithm samples the input data value, it does not integrate it. Thus total data value in the input image will not, in general, be conserved. However, an option is provided (see the " Control Flags" section below) which can produce approximate flux conservation by scaling the output values using the ratio of the output pixel size to the input pixel size. However, if accurate flux conservation is important to you, consder using the AST_REBIN$<$X$>$ or AST_REBINSEQ$<$X$>$ family of routines instead.

Output pixel coordinates are transformed into the coordinate system of the input grid using the inverse transformation of the Mapping which is supplied. This means that geometrical features in the input data are subjected to the Mapping’ s forward transformation as they are transferred from the input to the output grid (although the Mapping’ s forward transformation is not explicitly used).

In practice, transforming the coordinates of every pixel of a large data grid can be time-consuming, especially if the Mapping involves complicated functions, such as sky projections. To improve performance, it is therefore possible to approximate non-linear Mappings by a set of linear transformations which are applied piece-wise to separate sub-regions of the data. This approximation process is applied automatically by an adaptive algorithm, under control of an accuracy criterion which expresses the maximum tolerable geometrical distortion which may be introduced, as a fraction of a pixel.

This algorithm first attempts to approximate the Mapping with a linear transformation applied over the whole region of the output grid which is being used. If this proves to be insufficiently accurate, the output region is sub-divided into two along its largest dimension and the process is repeated within each of the resulting sub-regions. This process of sub-division continues until a sufficiently good linear approximation is found, or the region to which it is being applied becomes too small (in which case the original Mapping is used directly).

RESULT = AST_RESAMPLE$<$X$>$( THIS, NDIM_IN, LBND_IN, UBND_IN, IN, IN_VAR, INTERP, FINTERP, PARAMS, FLAGS, TOL, MAXPIX, BADVAL, NDIM_OUT, LBND_OUT, UBND_OUT, LBND, UBND, OUT, OUT_VAR, STATUS )

THIS = INTEGER (Given)

NDIM_IN = INTEGER (Given)

LBND_IN( NDIM_IN ) = INTEGER (Given)

UBND_IN( NDIM_IN ) = INTEGER (Given)

IN( $\ast$ ) = $<$Xtype$>$ (Given)

IN_VAR( $\ast$ ) = $<$Xtype$>$ (Given)

INTERP = INTEGER (Given)

FINTERP = SUBROUTINE (Given)

PARAMS( $\ast$ ) = DOUBLE PRECISION (Given)

FLAGS = INTEGER (Given)

TOL = DOUBLE PRECISION (Given)

MAXPIX = INTEGER (Given)

BADVAL = $<$Xtype$>$ (Given)

NDIM_OUT = INTEGER (Given)

LBND_OUT( NDIM_OUT ) = INTEGER (Given)

UBND_OUT( NDIM_OUT ) = INTEGER (Given)

LBND( NDIM_OUT ) = INTEGER (Given)

UBND( NDIM_OUT ) = INTEGER (Given)

OUT( $\ast$ ) = $<$Xtype$>$ (Returned)

OUT_VAR( $\ast$ ) = $<$Xtype$>$ (Returned)

STATUS = INTEGER (Given and Returned)

AST_RESAMPLE$<$X$>$ = INTEGER

• A value of zero will be returned if this function is invoked with the global error status set, or if it should fail for any reason.

To select the appropriate resampling function, you should replace $<$X$>$ in the generic function name AST_RESAMPLE$<$X$>$ with a 1- or 2-character data type code, so as to match the numerical type $<$Xtype$>$ of the data you are processing, as follows:

• D: DOUBLE PRECISION

• R: REAL

• I: INTEGER

• UI: INTEGER (treated as unsigned)

• S: INTEGER$\ast$2 (short integer)

• US: INTEGER$\ast$2 (short integer, treated as unsigned)

• B: BYTE (treated as signed)

• UB: BYTE (treated as unsigned)

For example, AST_RESAMPLED would be used to process DOUBLE PRECISION data, while AST_RESAMPLES would be used to process short integer data (stored in an INTEGER$\ast$2 array), etc.

For compatibility with other Starlink facilities, the codes W and UW are provided as synonyms for S and US respectively (but only in the Fortran interface to AST).

There is no such thing as a perfect sub-pixel interpolation scheme and, in practice, all resampling will result in some degradation of gridded data. A range of schemes is therefore provided, from which you can choose the one which best suits your needs.

In general, a balance must be struck between schemes which tend to degrade sharp features in the data by smoothing them, and those which attempt to preserve sharp features. The latter will often tend to introduce unwanted oscillations, typically visible as " ringing" around sharp features and edges, especially if the data are under-sampled (i.e. if the sharpest features are less than about two pixels across). In practice, a good interpolation scheme is likely to be a compromise and may exhibit some aspects of both these features.

For under-sampled data, some interpolation schemes may appear to preserve data resolution because they transform single input pixels into single output pixels, rather than spreading their data between several output pixels. While this may look better cosmetically, it can result in a geometrical shift of sharp features in the data. You should beware of this if you plan to use such features (e.g.) for image alignment.

The following are two easy-to-use sub-pixel interpolation schemes which are generally applicable. They are selected by supplying the appropriate value (defined in the AST_PAR include file) via the INTERP argument. In these cases, the FINTERP and PARAMS arguments are not used:

• AST__NEAREST: This is the simplest possible scheme, in which the value of the input pixel with the nearest centre to the interpolation point is used. This is very quick to execute and will preserve single-pixel features in the data, but may displace them by up to half their width along each dimension. It often gives a good cosmetic result, so is useful for quick-look processing, but is unsuitable if accurate geometrical transformation is required.

• AST__LINEAR: This is the default scheme, which uses linear interpolation between the nearest neighbouring pixels in the input grid (there are two neighbours in one dimension, four neighbours in two dimensions, eight in three dimensions, etc.). It is superior to the nearest-pixel scheme (above) in not displacing features in the data, yet it still executes fairly rapidly. It is generally a safe choice if you do not have any particular reason to favour another scheme, since it cannot introduce oscillations. However, it does introduce some spatial smoothing which varies according to the distance of the interpolation point from the neighbouring pixels. This can degrade the shape of sharp features in the data in a position-dependent way. It may also show in the output variance grid (if used) as a pattern of stripes or fringes.

An alternative set of interpolation schemes is based on forming the interpolated value from the weighted sum of a set of surrounding pixel values (not necessarily just the nearest neighbours). This approach has its origins in the theory of digital filtering, in which interpolated values are obtained by conceptually passing the sampled data (represented by a grid of delta functions) through a linear filter which implements a convolution. Because the convolution kernel is continuous, the convolution yields a continuous function which may then be evaluated at fractional pixel positions. The (possibly multi-dimensional) kernel is usually regarded as " separable" and formed from the product of a set of identical 1-dimensional kernel functions, evaluated along each dimension. Different interpolation schemes are then distinguished by the choice of this 1-dimensional interpolation kernel. The number of surrounding pixels which contribute to the result may also be varied.

From a practical standpoint, it is useful to divide the weighted sum of pixel values by the sum of the weights when determining the interpolated value. Strictly, this means that a true convolution is no longer being performed. However, the distinction is rarely important in practice because (for slightly subtle reasons) the sum of weights is always approximately constant for good interpolation kernels. The advantage of this technique, which is used here, is that it can easily accommodate missing data and tends to minimise unwanted oscillations at the edges of the data grid.

In the following schemes, which are based on a 1-dimensional interpolation kernel, the first element of the PARAMS array should be used to specify how many pixels are to contribute to the interpolated result on either side of the interpolation point in each dimension (the nearest integer value is used). Execution time increases rapidly with this number. Typically, a value of 2 is appropriate and the minimum value used will be 1 (i.e. two pixels altogether, one on either side of the interpolation point). A value of zero or less may be given for PARAMS(1) to indicate that a suitable number of pixels should be calculated automatically.

In each of these cases, the FINTERP argument is not used:

• AST__GAUSS: This scheme uses a kernel of the form exp(-k$\ast$x$\ast$x), with k a positive constant. The full-width at half-maximum (FWHM) is given by PARAMS(2) value, which should be at least 0.1 (in addition, setting PARAMS(1) to zero will select the number of contributing pixels so as to utilise the width of the kernel out to where the envelope declines to 1% of its maximum value). This kernel suppresses noise at the expense of smoothing the output array.

• AST__SINC: This scheme uses a sinc(pi$\ast$x) kernel, where x is the pixel offset from the interpolation point and sinc(z)=sin(z)/z. This sometimes features as an " optimal" interpolation kernel in books on image processing. Its supposed optimality depends on the assumption that the data are band-limited (i.e. have no spatial frequencies above a certain value) and are adequately sampled. In practice, astronomical data rarely meet these requirements. In addition, high spatial frequencies are often present due (e.g.) to image defects and cosmic ray events. Consequently, substantial ringing can be experienced with this kernel. The kernel also decays slowly with distance, so that many surrounding pixels are required, leading to poor performance. Abruptly truncating it, by using only a few neighbouring pixels, improves performance and may reduce ringing (if PARAMS(1) is set to zero, then only two pixels will be used on either side). However, a more gradual truncation, as implemented by other kernels, is generally to be preferred. This kernel is provided mainly so that you can convince yourself not to use it!

• AST__SINCSINC: This scheme uses an improved kernel, of the form sinc(pi$\ast$x).sinc(k$\ast$pi$\ast$x), with k a constant, out to the point where sinc(k$\ast$pi$\ast$x) goes to zero, and zero beyond. The second sinc() factor provides an " envelope" which gradually rolls off the normal sinc(pi$\ast$x) kernel at large offsets. The width of this envelope is specified by giving the number of pixels offset at which it goes to zero by means of the PARAMS(2) value, which should be at least 1.0 (in addition, setting PARAMS(1) to zero will select the number of contributing pixels so as to utilise the full width of the kernel, out to where it reaches zero). The case given by PARAMS(1)=2, PARAMS(2)=2 is typically a good choice and is sometimes known as the Lanczos kernel. This is a valuable general-purpose interpolation scheme, intermediate in its visual effect on images between the AST__NEAREST and AST__LINEAR schemes. Although the kernel is slightly oscillatory, ringing is adequately suppressed if the data are well sampled.

• AST__SINCCOS: This scheme uses a kernel of the form sinc(pi$\ast$x).cos(k$\ast$pi$\ast$x), with k a constant, out to the point where cos(k$\ast$pi$\ast$x) goes to zero, and zero beyond. As above, the cos() factor provides an envelope which gradually rolls off the sinc() kernel at large offsets. The width of this envelope is specified by giving the number of pixels offset at which it goes to zero by means of the PARAMS(2) value, which should be at least 1.0 (in addition, setting PARAMS(1) to zero will select the number of contributing pixels so as to utilise the full width of the kernel, out to where it reaches zero). This scheme gives similar results to the AST__SINCSINC scheme, which it resembles.

• AST__SINCGAUSS: This scheme uses a kernel of the form sinc(pi$\ast$x).exp(-k$\ast$x$\ast$x), with k a positive constant. Here, the sinc() kernel is rolled off using a Gaussian envelope which is specified by giving its full-width at half-maximum (FWHM) by means of the PARAMS(2) value, which should be at least 0.1 (in addition, setting PARAMS(1) to zero will select the number of contributing pixels so as to utilise the width of the kernel out to where the envelope declines to 1% of its maximum value). On astronomical images and spectra, good results are often obtained by approximately matching the FWHM of the envelope function, given by PARAMS(2), to the point spread function of the input data. However, there does not seem to be any theoretical reason for this.

• AST__SOMB: This scheme uses a somb(pi$\ast$x) kernel (a " sombrero" function), where x is the pixel offset from the interpolation point and somb(z)=2$\ast$J1(z)/z (J1 is a Bessel function of the first kind of order 1). It is similar to the AST__SINC kernel, and has the same parameter usage.

• AST__SOMBCOS: This scheme uses a kernel of the form somb(pi$\ast$x).cos(k$\ast$pi$\ast$x), with k a constant, out to the point where cos(k$\ast$pi$\ast$x) goes to zero, and zero beyond. It is similar to the AST__SINCCOS kernel, and has the same parameter usage.

In addition, the following schemes are provided which are not based on a 1-dimensional kernel:

• AST__BLOCKAVE: This scheme simply takes an average of all the pixels on the input grid in a cube centred on the interpolation point. The number of pixels in the cube is determined by the value of the first element of the PARAMS array, which gives the number of pixels in each dimension on either side of the central point. Hence a block of (2 $\ast$ PARAMS(1))$\ast$$\ast$NDIM_IN pixels in the input grid will be examined to determine the value of the output pixel. If the variance is not being used (USEVAR = .FALSE.) then all valid pixels in this cube will be averaged in to the result with equal weight. If variances are being used, then each input pixel will be weighted proportionally to the reciprocal of its variance; any pixel without a valid variance will be discarded. This scheme is suitable where the output grid is much coarser than the input grid; if the ratio of pixel sizes is R then a suitable value of PARAMS(1) may be R/2.

Finally, supplying the following values for INTERP allows you to implement your own sub-pixel interpolation scheme by means of your own routine. You should supply the name of this routine via the FINTERP argument:

• AST__UKERN1: In this scheme, you supply a routine to evaluate your own 1-dimensional interpolation kernel, which is then used to perform sub-pixel interpolation (as described above). The routine you supply should have the same interface as the fictitious AST_UKERN1 routine (q.v.). In addition, a value should be given via PARAMS(1) to specify the number of neighbouring pixels which are to contribute to each interpolated value (in the same way as for the pre-defined interpolation schemes described above). Other elements of the PARAMS array are available to pass values to your interpolation routine.

• AST__UINTERP: This is a completely general scheme, in which your interpolation routine has access to all of the input data. This allows you to implement any interpolation algorithm you choose, which could (for example) be non-linear, or adaptive. In this case, the AST_RESAMPLE$<$X$>$ functions play no role in the sub-pixel interpolation process and simply handle the geometrical transformation of coordinates and other housekeeping. The routine you supply should have the same interface as the fictitious AST_UINTERP routine (q.v.). In this case, the PARAMS argument is not used by AST_RESAMPLE$<$X$>$, but is available to pass values to your interpolation routine.

The following flags are defined in the AST_PAR include file and may be used to provide additional control over the resampling process. Having selected a set of flags, you should supply the sum of their values via the FLAGS argument:

• AST__NOBAD: Indicates that any output array elements for which no resampled value could be obtained should be left set to the value they had on entry to this function. If this flag is not supplied, such output array elements are set to the value supplied for argument BADVAL. Note, this flag cannot be used in conjunction with the AST__CONSERVEFLUX flag (an error will be reported if both flags are specified).

• AST__URESAMP1, 2, 3 & 4: A set of four flags which are reserved for your own use. They may be used to pass private information to any sub-pixel interpolation routine which you implement yourself. They are ignored by all the pre-defined interpolation schemes.

• AST__USEBAD: Indicates that there may be bad pixels in the input array(s) which must be recognised by comparing with the value given for BADVAL and propagated to the output array(s). If this flag is not set, all input values are treated literally and the BADVAL value is only used for flagging output array values.

• AST__USEVAR: Indicates that variance information should be processed in order to provide estimates of the statistical error associated with the resampled values. If this flag is not set, no variance processing will occur and the IN_VAR and OUT_VAR arrays will not be used. (Note that this flag is only available in the Fortran interface to AST.)

• AST__CONSERVEFLUX: Indicates that the output pixel values should be scaled in such a way as to preserve (approximately) the total data value in a feature on the sky. Without this flag, each output pixel value represents an instantaneous sample of the input data values at the corresponding input position. This is appropriate if the input data represents the spatial density of some quantity (e.g. surface brightness in Janskys per square arc-second) because the output pixel values will have the same normalisation and units as the input pixel values. However, if the input data values represent flux (or some other physical quantity) per pixel, then the AST__CONSERVEFLUX flag could be used. This causes each output pixel value to be scaled by the ratio of the output pixel size to the input pixel size.

This flag can only be used if the Mapping is successfully approximated by one or more linear transformations. Thus an error will be reported if it used when the TOL argument is set to zero (which stops the use of linear approximations), or if the Mapping is too non-linear to be approximated by a piece-wise linear transformation. The ratio of output to input pixel size is evaluated once for each panel of the piece-wise linear approximation to the Mapping, and is assumed to be constant for all output pixels in the panel. The scaling factors for adjacent panels will in general differ slightly, and so the joints between panels may be visible when viewing the output image at high contrast. If this is a problem, reduce the value of the TOL argument until the difference between adjacent panels is sufficiently small to be insignificant.

Note, this flag cannot be used in conjunction with the AST__NOBAD flag (an error will be reported if both flags are specified).

Unless the AST__NOBAD flag is specified, instances of missing data (bad pixels) in the output grid are identified by occurrences of the BADVAL value in the OUT array. These may be produced if any of the following happen:

• The input position (the transformed position of the output pixel’ s centre) lies outside the boundary of the grid of input pixels.

• The input position lies inside the boundary of a bad input pixel. In this context, an input pixel is considered bad if its data value is equal to BADVAL and the AST__USEBAD flag is set via the FLAGS argument. (Positions which have half-integral coordinate values, and therefore lie on a pixel boundary, are regarded as lying within the pixel with the larger, i.e. more positive, index.)

• The set of neighbouring input pixels (excluding those which are bad) is unsuitable for calculating an interpolated value. Whether this is true may depend on the sub-pixel interpolation scheme in use.

• The interpolated value lies outside the range which can be represented using the data type of the OUT array.

In addition, associated output variance estimates (if calculated) may be declared bad and flagged with the BADVAL value in the OUT_VAR array under any of the following circumstances:

• The associated resampled data value (in the OUT array) is bad.

• The set of neighbouring input pixels which contributed to the output data value do not all have valid variance estimates associated with them. In this context, an input variance estimate may be regarded as bad either because it has the value BADVAL (and the AST__USEBAD flag is set), or because it is negative.

• The set of neighbouring input pixels for which valid variance values are available is unsuitable for calculating an overall variance value. Whether this is true may depend on the sub-pixel interpolation scheme in use.

• The variance value lies outside the range which can be represented using the data type of the OUT_VAR array.

If the AST__NOBAD flag is specified via argument FLAGS, then output array elements that would otherwise be set to BADVAL are instead left holding the value they had on entry to this function. The number of such array elements is returned as the function value.

If the input or output grid is so large that an integer pixel index, (or a count of pixels) could exceed the largest value that can be represented by a 4-byte integer, then the alternative " 8-byte" interface for this function should be used. This alternative interface uses 8 byte integer arguments (instead of 4-byte) to hold pixel indices and pixel counts. Specifically, the arguments LBND_IN, UBND_IN, LBND_OUT, UBND_OUT, LBND, UBND are changed from type INTEGER to type INTEGER$\ast$8. The function return type is similarly changed to type INTEGER$\ast$8. The function name is changed by inserting the digit " 8" before the trailing data type code. Thus, AST_RESAMPLE$<$X$>$ becomes AST_RESAMPLE8$<$X$>$.