AST_RESAMPLE <X >

Resample a region of a data grid

Description:

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).

Invocation

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 )

Arguments

THIS = INTEGER (Given)
Pointer to a Mapping, whose inverse transformation will be used to transform the coordinates of pixels in the output grid into the coordinate system of the input grid. This yields the positions which are used to obtain resampled values by sub-pixel interpolation within the input grid.

The number of input coordinates used by this Mapping (as given by its Nin attribute) should match the number of input grid dimensions given by the value of NDIM_IN below. Similarly, the number of output coordinates (Nout attribute) should match the number of output grid dimensions given by NDIM_OUT.

NDIM_IN = INTEGER (Given)
The number of dimensions in the input grid. This should be at least one.
LBND_IN( NDIM_IN ) = INTEGER (Given)
An array containing the coordinates of the centre of the first pixel in the input grid along each dimension.
UBND_IN( NDIM_IN ) = INTEGER (Given)
An array containing the coordinates of the centre of the last pixel in the input grid along each dimension.

Note that LBND_IN and UBND_IN together define the shape and size of the input grid, its extent along a particular (J th) dimension being UBND_IN(J)-LBND_IN(J)+1. They also define the input grid s coordinate system, each pixel having unit extent along each dimension with integral coordinate values at its centre.

IN( ) = <Xtype > (Given)
An array, with one element for each pixel in the input grid, containing the input data to be resampled. The numerical type of this array should match the 1- or 2-character type code appended to the function name (e.g. if you are using AST_RESAMPLER, the type of each array element should be REAL).

The storage order of data within this array should be such that the index of the first grid dimension varies most rapidly and that of the final dimension least rapidly (i.e. normal Fortran array storage order).

IN_VAR( ) = <Xtype > (Given)
An optional second array with the same size and type as the IN array. If the AST__USEVAR flag is set via the FLAGS argument (below), this array should contain a set of non-negative values which represent estimates of the statistical variance associated with each element of the IN array. Estimates of the variance of the resampled output data will then be calculated.

If the AST__USEVAR flag is not set, no input variance estimates are required and this array will not be used. A dummy (e.g. one-element) array may then be supplied.

INTERP = INTEGER (Given)
This argument specifies the scheme to be used for sub-pixel interpolation within the input grid. It may be used to select from a set of pre-defined schemes by supplying one of the values described in the " Sub-Pixel Interpolation Schemes" section below. If a value of zero is supplied, then the default linear interpolation scheme is used (equivalent to supplying the value AST__LINEAR).

Alternatively, you may supply a value which indicates that you will provide your own routine to perform sub-pixel interpolation by means of the FINTERP argument. Again, see the " Sub-Pixel Interpolation Schemes" section below for details.

FINTERP = SUBROUTINE (Given)
If the value given for the INTERP argument indicates that you will provide your own routine for sub-pixel interpolation, then the name of that routine should be given here (the name should also appear in a Fortran EXTERNAL statement in the routine which invokes AST_RESAMPLE <X >). For details of the interface which the routine should have (several are possible, depending on the value of INTERP), see the " Sub-Pixel Interpolation Schemes" section below.

If the INTERP argument has any other value, corresponding to one of the pre-defined interpolation schemes, then this routine will not be used and you may supply the null routine AST_NULL here (note only one underscore). No EXTERNAL statement is required for this routine, so long as the AST_PAR include file has been used.

PARAMS( ) = DOUBLE PRECISION (Given)
An optional array which should contain any additional parameter values required by the sub-pixel interpolation scheme. If such parameters are required, this will be noted in the " Sub-Pixel Interpolation Schemes" section below (you may also use this array to pass values to your own interpolation routine).

If no additional parameters are required, this array is not used. A dummy (e.g. one-element) array may then be supplied.

FLAGS = INTEGER (Given)
The sum of a set of flag values which may be used to provide additional control over the resampling operation. See the " Control Flags" section below for a description of the options available. If no flag values are to be set, a value of zero should be given.
TOL = DOUBLE PRECISION (Given)
The maximum tolerable geometrical distortion which may be introduced as a result of approximating non-linear Mappings by a set of piece-wise linear transformations. This should be expressed as a displacement in pixels in the input grid s coordinate system.

If piece-wise linear approximation is not required, a value of zero may be given. This will ensure that the Mapping is used without any approximation, but may increase execution time.

MAXPIX = INTEGER (Given)
A value which specifies an initial scale size (in pixels) for the adaptive algorithm which approximates non-linear Mappings with piece-wise linear transformations. Normally, this should be a large value (larger than any dimension of the region of the output grid being used). In this case, a first attempt to approximate the Mapping by a linear transformation will be made over the entire output region.

If a smaller value is used, the output region will first be divided into sub-regions whose size does not exceed MAXPIX pixels in any dimension. Only at this point will attempts at approximation commence.

This value may occasionally be useful in preventing false convergence of the adaptive algorithm in cases where the Mapping appears approximately linear on large scales, but has irregularities (e.g. holes) on smaller scales. A value of, say, 50 to 100 pixels can also be employed as a safeguard in general-purpose software, since the effect on performance is minimal.

If too small a value is given, it will have the effect of inhibiting linear approximation altogether (equivalent to setting TOL to zero). Although this may degrade performance, accurate results will still be obtained.

BADVAL = <Xtype > (Given)
This argument should have the same type as the elements of the IN array. It specifies the value used to flag missing data (bad pixels) in the input and output arrays.

If the AST__USEBAD flag is set via the FLAGS argument, then this value is used to test for bad pixels in the IN (and IN_VAR) array(s).

Unless the AST__NOBAD flag is set via the FLAGS argument, this value is also used to flag any output elements in the OUT (and OUT_VAR) array(s) for which resampled values could not be obtained (see the " Propagation of Missing Data" section below for details of the circumstances under which this may occur). The AST_RESAMPLE <X > function return value indicates whether any such values have been produced. If the AST__NOBAD flag is set. then output array elements for which no resampled value could be obtained are left set to the value they had on entry to this function.

NDIM_OUT = INTEGER (Given)
The number of dimensions in the output grid. This should be at least one. It need not necessarily be equal to the number of dimensions in the input grid.
LBND_OUT( NDIM_OUT ) = INTEGER (Given)
An array containing the coordinates of the centre of the first pixel in the output grid along each dimension.
UBND_OUT( NDIM_OUT ) = INTEGER (Given)
An array containing the coordinates of the centre of the last pixel in the output grid along each dimension.

Note that LBND_OUT and UBND_OUT together define the shape, size and coordinate system of the output grid in the same way as LBND_IN and UBND_IN define the shape, size and coordinate system of the input grid.

LBND( NDIM_OUT ) = INTEGER (Given)
An array containing the coordinates of the first pixel in the region of the output grid for which a resampled value is to be calculated.
UBND( NDIM_OUT ) = INTEGER (Given)
An array containing the coordinates of the last pixel in the region of the output grid for which a resampled value is to be calculated.

Note that LBND and UBND together define the shape and position of a (hyper-)rectangular region of the output grid for which resampled values should be produced. This region should lie wholly within the extent of the output grid (as defined by the LBND_OUT and UBND_OUT arrays). Regions of the output grid lying outside this region will not be modified.

OUT( ) = <Xtype > (Returned)
An array, with one element for each pixel in the output grid, into which the resampled data values will be returned. The numerical type of this array should match that of the IN array, and the data storage order should be such that the index of the first grid dimension varies most rapidly and that of the final dimension least rapidly (i.e. normal Fortran array storage order).
OUT_VAR( ) = <Xtype > (Returned)
An optional array with the same type and size as the OUT array. If the AST__USEVAR flag is set via the FLAGS argument, this array will be used to return variance estimates for the resampled data values.

The output variance values will be calculated on the assumption that errors on the input data values are statistically independent and that their variance estimates may simply be summed (with appropriate weighting factors) when several input pixels contribute to an output data value. If this assumption is not valid, then the output error estimates may be biased. In addition, note that the statistical errors on neighbouring output data values (as well as the estimates of those errors) may often be correlated, even if the above assumption about the input data is correct, because of the sub-pixel interpolation schemes employed.

If the AST__USEVAR flag is not set, no output variance estimates will be calculated and this array will not be used. A dummy (e.g. one-element) array may then be supplied.

STATUS = INTEGER (Given and Returned)
The global status.

Returned Value

AST_RESAMPLE <X > = INTEGER
The number of output pixels for which no valid resampled value could be obtained. Thus, in the absence of any error, a returned value of zero indicates that all the required output pixels received valid resampled data values (and variances). See the BADVAL and FLAGS arguments.

Notes:

Data Type Codes

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:

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 INTEGER2 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).

Sub-Pixel Interpolation Schemes

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:

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:

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

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:

Control Flags

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:

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).

Propagation of Missing Data

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:

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:

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.

Handling of Huge Pixel Arrays

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 INTEGER8. The function return type is similarly changed to type INTEGER8. The function name is changed by inserting the digit " 8" before the trailing data type code. Thus, AST_RESAMPLE <X > becomes AST_RESAMPLE8 <X >.