This class grids and/or samples visibilities to/from UV space. More...

#include <wstackinggridder.h>

Public Types
enum	GridModeEnum { NearestNeighbourGridding, KaiserBesselKernel, RectangularKernel }
	Gridding modes that are supported for interpolating samples on the uv-grid. More...

Public Member Functions
	WStackingGridder (size_t width, size_t height, double pixelSizeX, double pixelSizeY, size_t fftThreadCount, ImageBufferAllocator *allocator, size_t kernelSize=7, size_t overSamplingFactor=63)
	Construct a new gridder with given settings. More...

	WStackingGridder (const WStackingGridder &)=delete

	~WStackingGridder ()
	De-allocate imagebuffers with the allocator and perform other clean-up. More...

void	AddData (const std::complex< float > *data, size_t dataDescId, double uInM, double vInM, double wInM)
	Grid an array of data values for inversion. More...

void	AddDataSample (std::complex< float > sample, double uInLambda, double vInLambda, double wInLambda)
	Grid a single visibility value for inversion. More...

ImageBufferAllocator *	Allocator () const
	The ImageBufferAllocator of this gridder. More...

void	FinalizeImage (double multiplicationFactor, bool correctFFTFactor)
	Finalize inversion once all passes are performed. More...

void	FinishInversionPass ()
	Finish an inversion gridding pass. More...

const std::complex< double > *	GetGriddedUVLayer (size_t layerIndex) const
	Retrieve a gridded uv layer. More...

void	GetGriddingCorrectionImage (double *image) const
	Make an image that contains the effect of the gridding kernel in image space. More...

void	GetKaiserBesselKernel (double *kernel, size_t n, bool multiplyWithSinc)
	Acquire a Kaiser-Bessel kernel. More...

enum GridModeEnum	GridMode () const
	Get the kernel function and its interpolation method used for gridding. More...

size_t	Height () const
	Get height of image as specified during construction. More...

double *	ImaginaryImage ()
	Get the imaginary part of a complex image after inversion. More...

void	InitializePrediction (const double *image)
	Initialize gridder for prediction and specify image to predict for. More...

void	InitializePrediction (const double real, const double imaginary)
	Complex alternative for InitializePrediction(const double*). More...

bool	IsComplex () const
	Whether the image produced by inversion or used by prediction is complex. More...

bool	IsInLayerRange (double wInLambda) const
	Determine whether specified w-value should be gridded in this pass. More...

bool	IsInLayerRange (double wStart, double wEnd) const
	Determine whether any samples within the specified w-value range should be gridded in this pass. More...

double	LayerToW (size_t layer) const
	Calculate the central w-value for the given w-layer index. More...

size_t	NFFTThreads () const
	Get the number of threads used when performing the FFTs. More...

size_t	NPasses () const
	Number of passes that are required when not all the w-layers fit in memory at once. More...

size_t	NWLayers () const
	Number of w-layers, as specified in PrepareWLayers(). More...

WStackingGridder &	operator= (const WStackingGridder &)=delete

double	PixelSizeX () const
	Get angular width of single pixel in radians. More...

double	PixelSizeY () const
	Get angular height of single pixel in radians. More...

void	PrepareBand (const MultiBandData &bandData)
	Initialize the inversion/prediction stage with a given band. More...

void	PrepareWLayers (size_t nWLayers, double maxMem, double minW, double maxW)
	Initialize the inversion/prediction stage with the given number of w-layers. More...

double *	RealImage ()
	Get the image result of inversion. More...

void	ReplaceImaginaryImageBuffer (double *newBuffer)
	Replace the current imaginary result of the imaging with a new image. More...

void	ReplaceRealImageBuffer (double *newBuffer)
	Replace the current result of the imaging with a new image. More...

void	SampleData (std::complex< float > *data, size_t dataDescId, double uInM, double vInM, double wInM)
	Predict the values for all channels in a band. More...

void	SampleDataSample (std::complex< double > &value, double uInLambda, double vInLambda, double wInLambda)
	Predict the value of a single visibility with double precision. More...

void	SampleDataSample (std::complex< float > &value, double uInLambda, double vInLambda, double wInLambda)
	Predict the value of a single visibility. More...

void	SetDenormalPhaseCentre (double dl, double dm)
	Setup the gridder for images with a shifted phase centre. More...

void	SetGridMode (enum GridModeEnum mode)
	Set the kernel function and its interpolation method used for gridding. More...

void	SetIsComplex (bool isComplex)
	Setup the gridder for real or complex images. More...

void	SetNFFTThreads (size_t nfftThreads)
	Set the number of threads used to perform the FFTs. More...

void	StartInversionPass (size_t passIndex)
	Initialize a new inversion gridding pass. More...

void	StartPredictionPass (size_t passIndex)
	Start a new prediction pass. More...

size_t	Width () const
	Get width of image as specified during construction. More...

size_t	WToLayer (double wInLambda) const
	Calculate on which layer this w-value belongs. More...

Detailed Description

This class grids and/or samples visibilities to/from UV space.

It also executes the FFT(s) required and performs the w-term correction. It does this with the 'w-stacking' method (see Offringa et al., 2014).

This class is quite generic in that it is agnostic to the underlying data format of the visibilities, and it is written to be of high performance. It is used by the WSMSGridder class to perform inversion and prediction. There, large imaging operations are split up in multiple 'passes' over the data if necessary, and in each pass a few w-layers are gridded/sampled. However, simpler applications might chose not to support this when very large images/large w-terms are not a concern, in which case things become easier.

This class can do both the calculation from visibilities to image, as well as going from (model) image to visibilities. The first is referred to as "Inversion" in this class, and the latter is called "Prediction".

The general call sequence for inversion (/imaging) is as follows:

Construct an instance with WStackingGridder()
Set settings if necessary (grid mode, denormal phase centre, complex...)
Call PrepareWLayers();
For each pass if multiple passes are necessary (or once otherwise) :
- Call StartInversionPass();
- Add all samples with AddDataSample();
- Call FinishInversionPass();
Finally, call FinalizeImage();
Now, RealImage() and optionally ImaginaryImage() will return the image(s).

Alternatively, AddData() can be used instead of AddDataSample, to grid several samples that only differ in frequency. To use AddData(), it is necessary to call PrepareBand() first.

For prediction, the sequence is similar:

Construct an instance with WStackingGridder()
Set settings if necessary
Call PrepareWLayers();
For each pass if multiple passes are necessary (or once otherwise) :
- Call InitializePrediction();
- Call StartPredictionPass();
- Get the predicted visibilities with SampleDataSample();

Similar to Inversion, all values for a band can be predicted at once by calling SampleData() instead of SampleDataSample() when the band has been set beforehand with PrepareBand().

Prediction does not require any finalisation calls.

Author: André Offringa

Date: 2013 (first version)

See also: WSClean: an implementation of a fast, generic wide-field imager for radio astronomy

Examples:: wspredictionexample.cpp.

Member Enumeration Documentation

enum WStackingGridder::GridModeEnum

Gridding modes that are supported for interpolating samples on the uv-grid.

Enumerator

NearestNeighbourGridding

Simple method that places/samples a visibility on the nearest uv-cell.

KaiserBesselKernel

Interpolate with a Kaiser-Bessel kernel.

This attenuates aliasing. The Kaiser-Bessel window is very similar to the prolate spheroidal kernel, which is considered the optimal gridding window. When this mode is selected, the kernel size and oversampling factor can be specified. This is the recommended and default mode.

RectangularKernel

Interpolate with a rectangular window.

This will give the sharpest transition at the edge of the image, so will maximally attenuate objects just past to the edge. However, objects further from the edge will not be as much attenuated compared to the KB window, which has much deeper sidelobes further out.

Constructor & Destructor Documentation

WStackingGridder::WStackingGridder	(	size_t	width,
		size_t	height,
		double	pixelSizeX,
		double	pixelSizeY,
		size_t	fftThreadCount,
		ImageBufferAllocator *	allocator,
		size_t	kernelSize = `7`,
		size_t	overSamplingFactor = `63`
	)

Construct a new gridder with given settings.

Currently, the width and height should be set equally.

Parameters

width	The width of the image in pixels
height	The height of the image in pixels.
pixelSizeX	The angular width of a pixel in radians.
pixelSizeY	The angular height of a pixel in radians.
fftThreadCount	The number of threads used for FFTs. See NFFTThreads() for more info.
allocator	An ImageBufferAllocator that is used for allocating images. See ImageBufferAllocator for the rational for using this.
kernelSize	The full width and height in pixels of the kernel. Should be odd. Large kernels cause smaller aliasing artefacts but are more expensive. The default of 7 is a reasonable compromise between effectiveness and speed.
overSamplingFactor	The number of different horizontal and vertical kernels that are precalculated at different rational positions. Larger is more accurate but requires more memory and becomes slower, probably mainly due to cache misses.

Todo:: Fix width/height requirement.

WStackingGridder::~WStackingGridder ( )

De-allocate imagebuffers with the allocator and perform other clean-up.

Member Function Documentation

void WStackingGridder::AddData	(	const std::complex< float > *	data,
		size_t	dataDescId,
		double	uInM,
		double	vInM,
		double	wInM
	)

Grid an array of data values for inversion.

The data values should have the same uvw-values in meters, i.e., they can only differ in frequency. The dataDescId specifies the frequencies of the array of data. This method requires that the channel frequencies have been specified beforehand, by calling PrepareBand().

This function will internally call AddDataSample() for each visibility, hence this function is just for convenience, but is not faster than individual calls to AddDataSample().

Parameters

data	Array of samples for different channels. The size of this array is given by the band referred to by dataDescId.
dataDescId	ID that specifies which band this data is for.
uInM	U value of UVW coordinate, in meters.
vInM	V value of UVW coordinate, in meters.
wInM	W value of UVW coordinate, in meters.

void WStackingGridder::AddDataSample	(	std::complex< float >	sample,
		double	uInLambda,
		double	vInLambda,
		double	wInLambda
	)

Grid a single visibility value for inversion.

This method does not require that the channel frequencies have been specified beforehand.

Parameters

sample	Visibility value.
uInLambda	U value of UVW coordinate, in number of wavelengths.
vInLambda	V value of UVW coordinate, in number of wavelengths.
wInLambda	W value of UVW coordinate, in number of wavelengths.

ImageBufferAllocator* WStackingGridder::Allocator ( ) const

inline

The ImageBufferAllocator of this gridder.

Returns: The image buffer allocator.

void WStackingGridder::FinalizeImage	(	double	multiplicationFactor,
		bool	correctFFTFactor
	)

Finalize inversion once all passes are performed.

Parameters

multiplicationFactor	Apply this factor to all pixels. This can be used to normalize the image for the weighting scheme.
correctFFTFactor	For normal imaging this should be false. It can be set to true when no correcting for weighting is performed, but the FFT factor should be taken out. WSClean uses this for certain values of "-visibility-weighting-mode".

void WStackingGridder::FinishInversionPass ( )

Finish an inversion gridding pass.

This will perform the Fourier transforms of the currently gridded w-layers, and add each gridded layer to the final image including w-term corrections. Therefore, it can take time.

See also: StartInversionPass().

const std::complex<double>* WStackingGridder::GetGriddedUVLayer ( size_t layerIndex ) const

inline

Retrieve a gridded uv layer.

This function can be called after StartInversionPass() was called, and before FinishInversionPass() is called.

Parameters

layerIndex Layer index of the grid, with zero being the first layer of the current pass.

Returns: The layer, with the currently gridded samples on it.

void WStackingGridder::GetGriddingCorrectionImage ( double * image ) const

Make an image that contains the effect of the gridding kernel in image space.

The gridder already corrects for the kernel internally, so it is not necessary to apply this image to the results, but users might be interested in the response function of the gridding kernel alone.

Parameters

image A width*height array that will be set to the kernel image.

void WStackingGridder::GetKaiserBesselKernel	(	double *	kernel,
		size_t	n,
		bool	multiplyWithSinc
	)

Acquire a Kaiser-Bessel kernel.

This is mostly a debugging/example function.

Parameters

kernel	Array of size `n`
n	Size of kernel
multiplyWithSinc	True to multiply the kernel with the sinc function, as is the case during gridding.

enum GridModeEnum WStackingGridder::GridMode ( ) const

inline

Get the kernel function and its interpolation method used for gridding.

Returns: The currently selected gridding mode.

size_t WStackingGridder::Height ( ) const

inline

Get height of image as specified during construction.

This is the full height, before any trimming has been applied.

Returns: Height in pixels of image being gridded.

double* WStackingGridder::ImaginaryImage ( )

inline

Get the imaginary part of a complex image after inversion.

Otherwise similar to RealImage().

void WStackingGridder::InitializePrediction ( const double * image )

inline

Initialize gridder for prediction and specify image to predict for.

PrepareWLayers() should be called before this method. Note that this method needs to be called with the same image before each pass, i.e., before each call to StartPredictionPass().

This method is used for prediction of non-complex (IsComplex()==false) images. Use InitializePrediction(const double*, const double*) for complex prediction – see SetIsComplex() for more info.

Parameters

image The model image that is to be predicted for. This is an array of width * height size, index by (x + width*y).

Examples:: wspredictionexample.cpp.

void WStackingGridder::InitializePrediction	(	const double *	real,
		const double *	imaginary
	)

inline

Complex alternative for InitializePrediction(const double*).

This method should be used for complex images – see SetIsComplex() for more info.

Parameters

real	Array of width * height giving the real (model) image values.
imaginary	Array of width * height giving the imaginary (model) image values.

bool WStackingGridder::IsComplex ( ) const

inline

Whether the image produced by inversion or used by prediction is complex.

In particular, cross-polarized images like XY and YX have complex values, because the UV grid is no longer Hermitian symmetric.

Returns: Whether the gridder is setup for complex images.

bool WStackingGridder::IsInLayerRange ( double wInLambda ) const

inline

Determine whether specified w-value should be gridded in this pass.

This method can only be called after calling StartInversionPass() or StartPredictionPass().

Parameters

wInLambda W-value in units of number of wavelengths.

Returns: true if the w-value is gridded on one of the w-layers that are processed in this pass.

bool WStackingGridder::IsInLayerRange	(	double	wStart,
		double	wEnd
	)		const

inline

Determine whether any samples within the specified w-value range should be gridded in this pass.

This can for example be used to optimize which samples to read from disc; when this function returns false for a given timestep/baseline, none of the samples are required in this pass. This method can only be called after calling StartInversionPass() or StartPredictionPass().

Parameters

wStart	W-value of range start in units of number of wavelengths.
wEnd	W-value of range end in units of number of wavelengths.

Returns: true if any of the w-values in the given range are processed in this pass.

double WStackingGridder::LayerToW ( size_t layer ) const

inline

Calculate the central w-value for the given w-layer index.

This is the inverse of WToLayer().

Parameters

layer W-layer index

Returns: Central w-value for specified w-layer in units of number of wavelengths.

size_t WStackingGridder::NFFTThreads ( ) const

inline

Get the number of threads used when performing the FFTs.

The w-layers are divided over the threads such that each thread calculates approximately the same number of w-layers. Therefore, optimally this value equals the number of (free) CPU cores available. If the number of w-layers is smaller than the number of threads, some threads will do nothing. Therefore, when not performing any w-term correction with NWLayers() == 1, this value has no effect.

When NWLayers() == 1 , it is still possible to perform multi-threading by setting up fftw to do so. The FFTWMultiThreadEnabler class can be used to setup fftw properly to do so.

Returns: The number of threads used for the FFTs.

size_t WStackingGridder::NPasses ( ) const

inline

Number of passes that are required when not all the w-layers fit in memory at once.

Valid once PrepareWLayers() has been called.

Returns: Number of passes.

Examples:: wspredictionexample.cpp.

size_t WStackingGridder::NWLayers ( ) const

inline

Number of w-layers, as specified in PrepareWLayers().

Returns: Number of w-layers.

double WStackingGridder::PixelSizeX ( ) const

inline

Get angular width of single pixel in radians.

The full image has a width of Width() * PixelSizeX(), assuming small angle approximation.

Returns: Width of a single pixel in radians.

double WStackingGridder::PixelSizeY ( ) const

inline

Get angular height of single pixel in radians.

The full image has a height of Height() * PixelSizeX(), assuming small angle approximation.

Returns: Height of a single pixel in radians.

void WStackingGridder::PrepareBand ( const MultiBandData & bandData )

inline

Initialize the inversion/prediction stage with a given band.

This is required for calling methods that take a dataDescId, specifically AddData() and SampleData(). This call can be avoided by using the AddDataSample() and SampleDataSample methods instead.

Parameters

bandData The information about the bands. A MultiBandData links a dataDescId to a set of contiguous frequencies. This corresponds with the DATA_DESC_ID field in meaurement sets.

void WStackingGridder::PrepareWLayers	(	size_t	nWLayers,
		double	maxMem,
		double	minW,
		double	maxW
	)

Initialize the inversion/prediction stage with the given number of w-layers.

The number of w-layers sets the accuracy of w-term correction and can normally be calculated from the image dimensions and maximum w-term.

After this call, NPasses() can be called to acquire the number of required passes.

Note that for non-complex images (the default), minW and maxW can be set to the minimum absolute w-value, which results in finer w-gridding and thus better w-term correction.

The minW and maxW values need to be calculated before gridding, and this normally requires an extra iteration through the w-values of the data before the data can be gridded.

For imaging without w-term correction, the number of layers can be set to one (note that one w-layer disables w-term correcting – do not set it to zero, as one uv-grid is always required). In that case, it is best to set minW to -maxW.

In cases where nwlayer > 1, the w-value w of all values should satisfy minW < abs(w) < maxW. When nWLayers == 1, this is not required.

Parameters

nWLayers	Number of uv grids at different w-values, should be >= 1.
maxMem	Allowed memory in bytes. The gridder will try to set the number of passes such that this value is not exceeded. Note that this is approximate.
minW	The smallest w-value to be inverted/predicted.
maxW	The largest w-value to be inverted/predicted.

Examples:: wspredictionexample.cpp.

double* WStackingGridder::RealImage ( )

inline

Get the image result of inversion.

This is an array of size width x height, and can be indexed with [x + y*width]. It is allowed to change this image, e.g. set the horizon to zero before saving to fits. This call is only valid once FinalizeImage() has been called.

If a complex image is produced, this image returns the real part. The imaginary part can be acquired with ImaginaryImage().

void WStackingGridder::ReplaceImaginaryImageBuffer ( double * newBuffer )

Replace the current imaginary result of the imaging with a new image.

Exactly like ReplaceRealImageBuffer(), except it replaces the imaginary part.

Parameters

newBuffer The new imaginary buffer part.

void WStackingGridder::ReplaceRealImageBuffer ( double * newBuffer )

Replace the current result of the imaging with a new image.

This is hardly ever really necessary but can save some memory and/or time when the WStackingGridder remains in memory with its results, and the results are passed on but need to be conditionally changed in some way (e.g. as done in the WSMSGridder class to pass a resized buffer to the WSClean class).

After calling this method, RealImage() will return newBuffer. The newBuffer will be cleaned up with the earlier provided ImageBufferAllocator. The old buffer can no longer be used.

This method replaces the real part of the imaging.

Parameters

newBuffer The new buffer that was allocated with the right ImageBufferAllocator.

See also: ReplaceImaginaryImageBuffer().

void WStackingGridder::SampleData	(	std::complex< float > *	data,
		size_t	dataDescId,
		double	uInM,
		double	vInM,
		double	wInM
	)

Predict the values for all channels in a band.

This is a convenience function that will call SampleDataSample() for all data values in the band specified by dataDescId.

If a particular value can not be predicted during this pass (because its w-value belongs to a w-layer of a different pass), its value will be set to NaN. A call to StartPredictionPass() should have been made before calling this method.

Parameters

data	Array of samples for different channels. The size of this array is given by the band referred to by dataDescId.
dataDescId	ID that specifies which band this data is for.
uInM	U value of UVW coordinate, in meters.
vInM	V value of UVW coordinate, in meters.
wInM	W value of UVW coordinate, in meters.

void WStackingGridder::SampleDataSample	(	std::complex< double > &	value,
		double	uInLambda,
		double	vInLambda,
		double	wInLambda
	)

Predict the value of a single visibility with double precision.

If the visibility can not be predicted, because its w-value belongs to a w-layer that is not processed during this pass, it will be given a value of NaN. A call to StartPredictionPass() should have been made before calling this method.

Parameters

value	Will be set to the predicted visibility value.
uInLambda	U value of UVW coordinate, in number of wavelengths.
vInLambda	V value of UVW coordinate, in number of wavelengths.
wInLambda	W value of UVW coordinate, in number of wavelengths.

Examples:: wspredictionexample.cpp.

void WStackingGridder::SampleDataSample	(	std::complex< float > &	value,
		double	uInLambda,
		double	vInLambda,
		double	wInLambda
	)

inline

Predict the value of a single visibility.

If the visibility can not be predicted, because its w-value belongs to a w-layer that is not processed during this pass, it will be given a value of NaN. A call to StartPredictionPass() should have been made before calling this method.

Parameters

value	Will be set to the predicted visibility value.
uInLambda	U value of UVW coordinate, in number of wavelengths.
vInLambda	V value of UVW coordinate, in number of wavelengths.
wInLambda	W value of UVW coordinate, in number of wavelengths.

void WStackingGridder::SetDenormalPhaseCentre	(	double	dl,
		double	dm
	)

inline

Setup the gridder for images with a shifted phase centre.

When dl or dm is non-zero, the centre of the involved images is shifted by the given amount. This allows inversion/prediction of images that are centred away from the phase centre, while the projection of the image is kept the same. The benefit of this can be that the phase centre can be set to zenith, thereby minimizing w-values, while the actually imaged field is somewhere else, thereby making it possible to handle non-zenith images with small image sizes and small w-values. A more extended explanation is given in the WSClean paper by Offringa et al. (2014).

Parameters

dl	shift of image centre in 'l' (East) direction
dm	shift of image centre in 'm' (North) direction

void WStackingGridder::SetGridMode ( enum GridModeEnum mode )

inline

Set the kernel function and its interpolation method used for gridding.

Parameters

mode	The new gridding mode.

void WStackingGridder::SetIsComplex ( bool isComplex )

inline

Setup the gridder for real or complex images.

See IsComplex().

Parameters

isComplex Whether the gridder should produce / predict from complex images.

void WStackingGridder::SetNFFTThreads ( size_t nfftThreads )

inline

Set the number of threads used to perform the FFTs.

See also: NFFTThreads() for an explanation.

Parameters

nfftThreads Number of threads used.

void WStackingGridder::StartInversionPass ( size_t passIndex )

Initialize a new inversion gridding pass.

PrepareWLayers() should have been called beforehand. Each call to StartInversionPass() should be followed by a call to FinishInversionPass().

Parameters

passIndex The zero-indexed index of this pass, 0 <= passIndex < NPasses().

void WStackingGridder::StartPredictionPass ( size_t passIndex )

Start a new prediction pass.

One of the InitializePrediction() methods should be called before each call to this method. This method will perform the fast Fourier transforms necessary for this pass, so it can take time.

Parameters

passIndex Zero-indexed index of prediction pass, 0 <= passIndex < NPasses().