Gridder utility functions

These functions are mainly internal utilities for the w-towers sub-grid gridder and degridder. The most useful ones are likely to be ska_sdp_func.grid_data.determine_w_step() and ska_sdp_func.grid_data.determine_max_w_tower_height(), which can be used to evaluate those parameters which the gridder requires.

All the Python functions call the underlying C versions, which should be faster than pure Python implementations.

C/C++

struct sdp_GridProjection

Represents the transformation of pixel coordinates into the original (l,m,n) coordinate system as aligned with the (u,v,w) axes associated with the incoming visibilities.

Our image coordinate system @(l’,m’)@ is defined as follows, with @px_l@ the coordinate of the first dimension and @px_m@ the coordinate of the second dimension:

\[\begin{split} \begin{align*} px_l &= size_l \left(\frac{l'}{\theta_l} + \frac12\right) \\ px_m &= size_m \left(\frac{m'}{\theta_m} + \frac12\right) \end{align*} \end{split}\]

This will map as follows to the original (l,m,n) sky coordinates relative to the visibility phase centre:

\[\begin{split} \begin{pmatrix} l' \\ m' \\ n' \end{pmatrix} = \begin{pmatrix} 1 & 0 & shear_u \\ 0 & 1 & shear_v \\ rectr_l & rectr_m & 1 \end{pmatrix} \begin{pmatrix} l \\ m \\ n \end{pmatrix} - \begin{pmatrix} shift_l \\ shift_m \\ shift_n \end{pmatrix} \end{split}\]

Where the following can be set by the user:

shift_l/shift_m are offsets to shift the centre of the image. Note that these get applied after the linear coordinate transformation (so they can be used to precisely line up sub-images even with a shear transformation applied).
shear_u/shear_v are a shear transformation that can be used to minimise $|w|$ in the @(u,v,w)@ layout (see below)

The remaining transformation parameters are set automatically, as they only affect $n’$ (which cannot be observed externally):

shift_n is determined in a way that $|n’|$ is minimised
rectr_l/rectrl_m are used internally to ce-centre visibilities on their main contribution to the grid. This maximises locality and minimises the $n’$ range we need to cover.

The corresponding UVW transformation that is applied interally is as follows (inverse transpose of above matrix):

\[\begin{split} \begin{pmatrix} u' \\ v' \\ w' \end{pmatrix} = \frac 1{1-shear_u rectr_l - shear_v rectr_m} \begin{pmatrix} 1-shear_v*rectr_m & shear_v*rectr_l & -rectr_l \\ shear_u*rectr_m & 1-shear_u*rectr_l & -rectr_m \\ -shear_u & -shear_v & 1 \end{pmatrix} \begin{pmatrix} u \\ v \\ w \end{pmatrix} \end{split}\]

Note that @w’ = w - shear_u u - shear_v v@, which means that if @shear_u@ and @shear_v@ are chosen suitably, we can minimise $w’$ (e.g. by linear regression over the UVW layout). This can greatly reduce the complexity of gridding, at the price of a - typically small - image warp.

Public Members

int64_t size_l: Size of image / grid in pixels along l axis.

int64_t size_m: Size of image / grid in pixels along m axis.

double theta_l: Size of image in direction cosines on l axis / resolution of grid on u axis

double theta_m: Size of image in direction cosines on m axis / resolution of grid on v axis

double shift_l: Shift to apply to phase (and therefore image) centre on l axis.

double shift_m: Shift to apply to phase (and therefore image) centre on m axis.

double shear_u: Shear to apply to visibility w positions depending on u axis (warping the image on l axis depending on n)

double shear_v: Shear to apply to visibility w positions depending on v axis (warping the image on m axis depending on n)

double w_step: Step length for w-towers. The resolution of our @(u’,v’,w’)@ grid will therefore be @(1/theta_l,1/theta_m,w_step).

double shift_n: Shift to apply to phase centre on n axis. In contrast to shift_l and shift_m, this shift is cancelled out inside the gridder. Or put another way: While this parameter indeed shifts the internal coordinate system as described, the only change in external behaviour will be performance and numerical accuracy. Check determine_n_shift for how to set this parameter optimally.

double rectr_l

Recentering factors. These are chosen such that partial derivates of the $n’$ coordinates corresponding to the sky sphere (i.e. \sqrt{1-l^2-m^2}@ after transformation) at @l’=m’=0@ are zero. This has a number of advantages - it means the maximum of

’@ is always in the centre of the image (which minimises the

’@ range to cover), and even more critically maximises the symmetry of

’@ and therefore the one of the @w-kernel, “recentering” it. In practical terms, this both maximises w-tower step length and total height, significantly boosting efficency when shifting the image.

void sdp_gridder_clamp_channels_single(const sdp_Mem *uvws, const int dim, const double freq0_hz, const double dfreq_hz, const sdp_Mem *start_ch_in, const sdp_Mem *end_ch_in, const double min_u, const double max_u, sdp_Mem *start_ch_out, sdp_Mem *end_ch_out, int64_t start_row, int64_t end_row, sdp_Error *status)

Clamp channels for a single dimension of an array of uvw coordinates.

Restricts a channel range such that all visibilities lie in the given range in u or v or w.

Parameters:

uvws – float[uvw_count, 3] UVW coordinates of visibilities (in m).
dim – Dimension index (0, 1 or 2) of uvws to check.
freq0_hz – Frequency of first channel, in Hz.
dfreq_hz – Channel width, in Hz.
start_ch_in – Input channel range to clamp (excluding end).
end_ch_in – Input channel range to clamp (excluding end).
min_u – Minimum value for u or v or w (inclusive).
max_u – Maximum value for u or v or w (exclusive).
start_ch_out – Clamped start channel.
end_ch_out – Clamped end channel (excluding end).
start_row – Row (uvw index) at which to start processing data.
end_row – Row (uvw index) at which to stop processing data (exclusive).
status – Error status.

void sdp_gridder_clamp_channels_uv(const sdp_Mem *uvws, const double freq0_hz, const double dfreq_hz, const sdp_Mem *start_ch_in, const sdp_Mem *end_ch_in, const double min_u, const double max_u, const double min_v, const double max_v, sdp_Mem *start_ch_out, sdp_Mem *end_ch_out, int64_t start_row, int64_t end_row, sdp_Error *status)

Clamp channels for (u,v) in an array of uvw coordinates.

Restricts a channel range such that all visibilities lie in the given range in u and v.

Parameters:

uvws – float[uvw_count, 3] UVW coordinates of visibilities (in m).
freq0_hz – Frequency of first channel, in Hz.
dfreq_hz – Channel width, in Hz.
start_ch_in – Input channel range to clamp (excluding end).
end_ch_in – Input channel range to clamp (excluding end).
min_u – Minimum value for u (inclusive).
max_u – Maximum value for u (exclusive).
min_v – Minimum value for v (inclusive).
max_v – Maximum value for v (exclusive).
start_ch_out – Clamped start channel.
end_ch_out – Clamped end channel (excluding end).
start_row – Row (uvw index) at which to start processing data.
end_row – Row (uvw index) at which to stop processing data (exclusive).
status – Error status.

Warning

doxygenfunction: Cannot find function “sdp_gridder_determine_w_step” in doxygen xml output for project “ska-sdp-func” from directory: ../doxygen/xml

int sdp_gridder_determine_max_w_tower_height(sdp_GridProjection proj, int subgrid_size, double w_max, int support, int oversampling, int w_support, int w_oversampling, double fov_l, double fov_m, double subgrid_frac, int num_samples, double target_err, sdp_Error *status)

Find the maximum w-tower height for a given configuration.

This function is ported from the implementation provided in Peter Wortmann’s subgrid_imaging Jupyter notebook.

Note that the image size supplied does not need to be the same as the final image - and for performance reasons you might want to use a smaller image size. Smaller images generally yield somewhat higher errors (the errors “wrap around”). Using about ~twice the subgrid size as the image size generally means that you will be underestimating the viable w-tower height by only a few percents.

Parameters:

proj – Image grid projection to assume
subgrid_size – Sub-grid size in pixels.
w_max – Maximum height to consider for w-towers.
support – Kernel support size in (u, v).
oversampling – Oversampling factor for uv-kernel.
w_support – Support size in w.
w_oversampling – Oversampling factor for w-kernel.
fov_l – Total image size in direction cosines along l-axis.
fov_m – Total image size in direction cosines along m-axis.
subgrid_frac – Fraction of subgrid size that should be used.
num_samples – Number of samples in u and v (if 0, defaults to 3).
target_err – Optional target error (if 0, aim for twice as much error as gridder has in image centre).
status – Error status.

void sdp_gridder_make_kernel(const sdp_Mem *window, sdp_Mem *kernel, sdp_Error *status)

Convert image-space window function to oversampled kernel.

This uses a DFT to do the transformation to Fourier space. The supplied input window function is 1-dimensional, with a shape of (support), and the output kernel is 2-dimensional, with a shape of (oversample + 1, support), and it must be sized appropriately on entry.

Parameters:

window – Image-space window function in 1D.
kernel – Oversampled output kernel in Fourier space.
status – Error status.

void sdp_gridder_make_pswf_kernel(int support, sdp_Mem *kernel, sdp_Error *status)

Generate an oversampled kernel using PSWF window function.

This uses a DFT to do the transformation to Fourier space. The output kernel is 2-dimensional, with a shape of (oversample + 1, vr_size), and it must be sized appropriately on entry.

Parameters:

support – Support size for kernel, usually the same as vr_size.
kernel – Oversampled output kernel in Fourier space.
status – Error status.

void sdp_gridder_make_w_pattern(sdp_GridProjection proj, sdp_Mem *w_pattern, sdp_Error *status)

Generate w-pattern.

This is the iDFT of a single visibility at (0, 0, w).

Parameters:

subgrid_size – Subgrid size.
theta – Total image size in direction cosines on l axis.
theta – Total image size in direction cosines on m axis.
shear_u – Shear parameter in u.
shear_v – Shear parameter in v.
w_step – Distance between w-planes.
w_pattern – Complex w-pattern, dimensions (subgrid_size, subgrid_size).
status – Error status.

double sdp_gridder_rms_diff(const sdp_Mem *a, const sdp_Mem *b, sdp_Error *status)

Returns the RMS of the difference between two 2D arrays: rms(a - b).

The two arrays must be 2D and have the same shape.

Parameters:

a – The first input array.
b – The second input array.
status – Error status.

void sdp_gridder_subgrid_add(sdp_Mem *grid, int offset_u, int offset_v, const sdp_Mem *subgrid, double factor, sdp_Error *status)

Add the supplied sub-grid to the grid.

Parameters:

grid – Output grid.
offset_u – Offset in u.
offset_v – Offset in v.
subgrid – Input sub-grid.
factor – Factor by which to multiply elements of sub-grid before adding.
status – Error status.

void sdp_gridder_subgrid_cut_out(const sdp_Mem *grid, int offset_u, int offset_v, sdp_Mem *subgrid, sdp_Error *status)

Cut out a sub-grid from the supplied grid.

Parameters:

grid – Input grid.
offset_u – Offset in u.
offset_v – Offset in v.
subgrid – Output sub-grid.
status – Error status.

void sdp_gridder_uvw_bounds_all(const sdp_Mem *uvws, double freq0_hz, double dfreq_hz, const sdp_Mem *start_chs, const sdp_Mem *end_chs, double uvw_min[3], double uvw_max[3], sdp_Error *status)

Determine (scaled) min and max values in uvw coordinates.

Parameters:

uvws – float[uvw_count, 3] UVW coordinates of visibilities (in m).
freq0_hz – Frequency of first channel (Hz).
dfreq_hz – Channel separation (Hz).
start_chs – int[uvw_count] First channel to degrid for every uvw.
end_chs – int[uvw_count] Channel at which to stop degridding for every uvw.
uvw_min – Output 3-element array containing minimum (u,v,w) values.
uvw_max – Output 3-element array containing maximum (u,v,w) values.
status – Error status.

Python

class ska_sdp_func.grid_data.GridProjection

Wrapper for GridProjection structure

Represents the transformation of pixel coordinates into the original (l,m,n) coordinate system as aligned with the (u,v,w) axes associated with the incoming visibilities. We assume that the following transformation has been applied to the original (l,m,n) coordinates:

(l')   (1  0  shear_u) (l)   (shift_l)
(m') = (0  1  shear_v) (m) - (shift_m)
(n')   (0  0  1      ) (n)   (shift_n)

Corresponding to a phase rotation of the visibilities for the visibilities (to implement the shift of the centre), and an UVW transformation of:

(u')    (1         0         0) (u)
(v') =  (0         1         0) (v)
(w')    (-shear_u  -shear_v  1) (w)

to implement the shear. Then image pixel coordinates follow as:

px_l = size_l (l' / theta_l + 0.5)
px_m = size_m (m' / theta_m + 0.5)

We are clearly typically interested in the opposite transformation, which is not quite obvious (see px_to_l, px_to_m and lm_to_n)

Parameters:

size_l – Pixel size of image / grid in pixels on l axis
size_m – Pixel size of image / grid in pixels on m axis
theta_l – Size of image in direction cosines / resolution of grid on l axis
theta_m – Size of image in direction cosines / resolution of grid on m axis
shift_l – Shift to apply to phase (and therefore image) centre on l axis
shift_m – Shift to apply to phase (and therefore image) centre on m axis
shear_u – Shear to apply to visibilities depending on u coordinate (warping the image)
shear_v – Shear to apply to visibilities depending on v coordinate (warping the image)

ska_sdp_func.grid_data.clamp_channels_single(uvws: ndarray, dim: int, freq0_hz: float, dfreq_hz: float, start_ch_in: ndarray, end_ch_in: ndarray, min_u: float, max_u: float, start_ch_out: ndarray, end_ch_out: ndarray, start_row: int = -1, end_row: int = -1)

Clamp channels for a single dimension of an array of uvw coordinates.

Restricts a channel range such that all visibilities lie in the given range in u or v or w.

Parameters:

uvws – float[uvw_count, 3] UVW coordinates of visibilities (in m).
dim – Dimension index (0, 1 or 2) of uvws to check.
freq0_hz – Frequency of first channel, in Hz.
dfreq_hz – Channel width, in Hz.
start_ch_in – Input channel range to clamp (excluding end).
end_ch_in – Input channel range to clamp (excluding end).
min_u – Minimum value for u or v or w (inclusive).
max_u – Maximum value for u or v or w (exclusive).
start_ch_out – Clamped start channel.
end_ch_out – Clamped end channel (excluding end).
start_row – Row (uvw index) at which to start processing data.
end_row – Row (uvw index) at which to stop processing data (exclusive).

ska_sdp_func.grid_data.clamp_channels_uv(uvws: ndarray, freq0_hz: float, dfreq_hz: float, start_ch_in: ndarray, end_ch_in: ndarray, min_u: float, max_u: float, min_v: float, max_v: float, start_ch_out: ndarray, end_ch_out: ndarray, start_row: int = -1, end_row: int = -1)

Clamp channels for (u,v) in an array of uvw coordinates.

Restricts a channel range such that all visibilities lie in the given range in u and v.

Parameters:

uvws – float[uvw_count, 3] UVW coordinates of visibilities (in m).
freq0_hz – Frequency of first channel, in Hz.
dfreq_hz – Channel width, in Hz.
start_ch_in – Input channel range to clamp (excluding end).
end_ch_in – Input channel range to clamp (excluding end).
min_u – Minimum value for u (inclusive).
max_u – Maximum value for u (exclusive).
min_v – Minimum value for v (inclusive).
max_v – Maximum value for v (exclusive).
start_ch_out – Clamped start channel.
end_ch_out – Clamped end channel (excluding end).
start_row – Row (uvw index) at which to start processing data.
end_row – Row (uvw index) at which to stop processing data (exclusive).

ska_sdp_func.grid_data.determine_max_w_tower_height(proj: GridProjection, subgrid_size: int, fov_l: float, fov_m: float, support: int, oversampling: int, w_support: int, w_oversampling: int, subgrid_frac: float = 0.6666666666666666, num_samples: int = 3, target_err: float | None = None, w_max: float = 0.0) → float

Find maximum w-tower height of a given configuration by trial-and-error.

This is the same as find_max_w_tower_height(), but without needing to create and supply a gridder kernel.

Parameters:

proj – Image grid projection. If image_size is 0, it will be set to twice the subgrid size.
subgrid_size – Sub-grid size in pixels.
theta_l – Total image size on l axis in direction cosines.
theta_m – Total image size on m axis in direction cosines.
fov – Un-padded image field of view on l axis, in direction cosines.
fov – Un-padded image field of view on m axis, in direction cosines.
support – Kernel support size in (u, v).
oversampling – Oversampling factor for uv-kernel.
w_support – Support size in w.
w_oversampling – Oversampling factor for w-kernel.
image_size – Size of test image, in pixels. If not specified, defaults to twice the subgrid size.
shear_u – Shear parameter in u (use zero for no shear).
shear_v – Shear parameter in v (use zero for no shear).
subgrid_frac – Fraction of sub-grid to use.
num_samples – Number of sample points to test in u and v directions.
target_err – Target error to use. If None, it is determined automatically.
w_max – Maximum w to consider. Providing this sort of hint can speed finding the right height up considerably.

Returns:

The maximum w-tower height in units of w_step.

ska_sdp_func.grid_data.find_max_w_tower_height(grid_kernel: GridderWtowerUVW, fov_l: float, fov_m: float, subgrid_frac: float = 0.6666666666666666, num_samples: int = 3, target_err: float | None = None, w_max: float = 0.0)

Find maximum w-tower height of a given configuration by trial-and-error.

This matches the interface with the function in the notebook.

Parameters:

grid_kernel – Gridder kernel to use for the evaluation.
fov_l – Un-padded image field of view on l axis, in direction cosines.
fov_m – Un-padded image field of view on m axis, in direction cosines.
subgrid_frac – Fraction of sub-grid to use.
num_samples – Number of sample points to test in u and v directions.
target_err – Target error to use. If None, it is determined automatically.
w_max – Maximum w to consider. Providing this sort of hint can speed finding the right height up considerably.

Returns:

The maximum w-tower height in units of grid_kernel.w_step.

ska_sdp_func.grid_data.make_kernel(window: ndarray, kernel: ndarray)

Convert image-space window function to oversampled kernel.

This uses a DFT to do the transformation to Fourier space. The supplied input window function is 1-dimensional, with a shape of (support), and the output kernel is 2-dimensional, with a shape of (oversample + 1, support), and it must be sized appropriately on entry.

Parameters:

window – Image-space window function in 1D.
kernel – Oversampled output kernel in Fourier space.

ska_sdp_func.grid_data.make_pswf_kernel(support: int, kernel: ndarray)

Generate an oversampled kernel using PSWF window function.

This uses a DFT to do the transformation to Fourier space. The output kernel is 2-dimensional, with a shape of (oversample + 1, support), and it must be sized appropriately on entry.

Parameters:

support – Kernel support size.
kernel – Oversampled output kernel in Fourier space.

ska_sdp_func.grid_data.make_w_pattern(proj: GridProjection) → ndarray

Generate w-pattern.

This is the iDFT of a single visibility at (0, 0, w).

Parameters:

proj – Image grid projection to assume.
w_step – Distance between w-planes.

Returns:

Complex w-pattern output, dimensions (size_l, size_m).

ska_sdp_func.grid_data.rms_diff(array_a: ndarray, array_b: ndarray)

Returns the RMS of the difference between two 2D arrays: rms(a - b).

The two arrays must be 2D and have the same shape.

Parameters:

array_a – The first input array.
array_b – The second input array.

ska_sdp_func.grid_data.subgrid_add(grid: ndarray, offset_u: int, offset_v: int, subgrid: ndarray, factor: float = 1.0)

Add the supplied sub-grid to the grid.

Parameters:

grid – Output grid.
offset_u – Offset in u.
offset_v – Offset in v.
subgrid – Input sub-grid.
factor – Factor by which to multiply elements of sub-grid before adding.

ska_sdp_func.grid_data.subgrid_cut_out(grid: ndarray, offset_u: int, offset_v: int, subgrid: ndarray)

Cut out a sub-grid from the supplied grid.

Parameters:

grid – Input grid.
offset_u – Offset in u.
offset_v – Offset in v.
subgrid – Output sub-grid.

ska_sdp_func.grid_data.uvw_bounds_all(uvws: ndarray, freq0_hz: float, dfreq_hz: float, start_ch: ndarray, end_ch: ndarray)

Determine (scaled) min and max values in uvw coordinates.

Parameters:

uvws – float[uvw_count, 3] UVW coordinates of visibilities (in m).
freq0_hz – Frequency of first channel, in Hz.
dfreq_hz – Channel width, in Hz.
start_ch – First channel to degrid for every uvw.
end_ch – Channel at which to stop degridding for every uvw.

Returns:

A tuple of two 3-element lists containing minimum and maximum (u,v,w)-values: (uvw_min, uvw_max)