The dimension of the convolution done is determined by the dimension of
the point spread function (psf), so for example, if the psf is a Vector,
one dimensional convolution will be done. The dimension of the model
that is to be convolved must be at least the same as the point
spread function, but it can be larger. If it is then the convolution will
be repeated for each row or plane of the model.
Repeated convolution can only be done along the fastest moving axes of
the supplied image. For example, if a one dimensional psf is used (so
that one dimensional convolution is being done), and a cube of data is
supplied then the convolution will be repeated for each row in the
cube. It is not currently possible to have this class do repeated one
dimensional convolution along all the columns or along the z
axis. To do this you need to use an iterator external to the class to
successively feed in the appropriate slices of your Array.
The difference between linear and circular convolution can best be
explained with a one dimensional example.
Suppose the psf and data to be convolved are:
The "centre" of the convolution is at the point (NX/2, NY/2) (assuming a
2 dimensional psf) where the first point in the psf is at (0,0) and the
last is at (NX-1, NY-1). This means that a psf that is all zero except
for 1 at the "centre" pixel will when convolved with any model leave the
model unchanged.
The convolution is done in the Fourier domain and the transform of the
psf (the transfer function) is cached by this class. If the cached
transfer function is the wrong size for a given model it will be
automatically be recomputed to the right size (this will involve two
FFT's)
Each convolution requires two Fourier transforms which dominate the
computational load. Hence the computational expense is
n Log(n) for 1 dimensional and
n^2 Log(n) for 2 dimensional convolutions.
The size of the convolved result is always the same as the input model
unless linear convolution is done with the fullSize option set to True.
In this case the result will be larger than the model and include the
full linear convolution (resultSize = psfSize+modelSize-1), rather than
the central portion.
If the convolver is constructed with an expected model size (as in the
example below) then the cached transfer function will be computed to a
size appropriate for linear convolution of models of that size. If no
model size is given then the cached transfer function will be computed
with a size appropriate for circular convolution. These guidelines also
apply when using the setPsf functions.
For linear convolution the psf can be larger, the same size or smaller
than the model but for circular convolution the psf must be smaller or the
same size.
The size of the cached transfer function (and also the length of the
FFT's calculated) depends on the sizes of the psf and the model, as well
as whether you are doing linear or circular convolution and the fullSize
option. It is always advantageous to use the smallest possible psf
(ie. do not pad the psf prior to supplying it to this class). Be aware
that using odd length images will lead to this class doing odd length
FFT's, which are less computationally effecient (particularly is the
length of the transform is a prime number) in general than even length
transforms.
There are only two valid template types namely,
When this class is constructed you may choose to have the psf
explicitly stored by the class (by setting cachePsf=True). This will
allow speedy access to the psf when using the getPsf function. However
the getPsf function can still be called even if the psf has not been
cached. Then the psf will be computed by FFT'ing the transfer
function, and the psf will also then be cached (unless
cachePsf=Flase). Cacheing the psf is also a good idea if you will be
switching between different sized transfer functions (eg. mixing
linear and circular convolution) as it will save one of the two
FFT's. Note that even though the psf is returned as a const Array, it
is possible to inadvertently modify it using the array copy constructor
as this uses reference symantics. Modifying the psf is NOT
recommended. eg.
The copy constructor and the assignment operator make copies (and not
references) of all the internal data arrays, as this object could get
really screwed up if the private data was silently messed with.
The destructor does nothing!
Perform linear convolution of the model with the previously
specified psf. Return the answer in result. Set fullSize to True if you
want the full convolution, rather than the central portion (the same
size as the model) returned.
Perform circular convolution of the model with the previously
specified psf. Return the answer in result.
Set the transfer function for future convolutions to psf.
Assume circular convolution will be done
Set to use convolution with lesser flips
Prerequisite
The mathematical concept of convolution
Etymology
The convolver class performs convolution!
Synopsis
This class will perform linear or circular convolution on arrays.
This class strips all degenerate axes when determining the dimensionality
of the psf or model. So a psf with shapes of [1,1,16] or [16,1,1] is
treated the same as a Vector of length 16, and will result in one
dimensional convolutions along the first non-degenerate axis of the
supplied model.
psf = [0 .5 1 .1]; data = [1 0 0 0 0 0]
then their linear and circular convolutions are:
circular convolution = [1 .1 0 0 0 .5]
linear convolution = [1 .1 0 0 0 0] (fullSize == False)
linear convolution = [0 .5 1 .1 0 0 0 0 0] (fullSize == True)
The circular convolution "wraps around" whereas the linear one does not.
Usage of the fullSize option is explained below. As can be seen from the
above example this class does not normalise the convolved result by any
factor that depends on the psf, so if the "beam area" is not unity the
flux scales will vary.
If you are intending to do 'fullsize' linear convolutions
you should also set the fullsize option to True as the cached transfer
function is a different size for fullsize linear convolutions.
and the user may prefer to use the following typedef's:
FloatConvolver (= Convolver<Float>) or
DoubleConvolver (= Convolver<Double>)
rather than explicitly specifying the template arguements.
The typedefs need to be redeclared when using the gnu compiler making
them essentially useless.
DoubleConvolver conv();
{
Matrix<Double> psf(20,20);
conv.setPsf(psf);
}
Matrix<Double> convPsf = conv.getPsf(); // Get the psf used by the convolver
convPsf(0,0) = -100; // And modify it. This modifies
// This internal psf used by the
// convolver also! (unless it is
// not caching the psf)
Example
Calculate the convolution of two Matrices (psf and model);
Matrix<Float> psf(4,4), model(12,12);
...put meaningful values into the above two matrices...
FloatConvolver conv(psf, model.shape());
conv.linearConv(result, model); // result = Convolution(psf, model)
Motivation
I needed to do linear convolution to write a clean algorithm. It
blossomed into this class.
Thrown Exceptions
To Do
Member Description
Convolver()
When using the default constructor the psf MUST be specified using the
setPsf function prior to doing any convolution.
Convolver(const Array<FType>& psf, Bool cachePsf=False)
Create the cached Transfer function assuming that circular convolution
will be done
Convolver(const Array<FType>& psf, const IPosition& imageSize, Bool fullSize=False, Bool cachePsf=False)
Create the cached Transfer function assuming that linear convolution
with an array of size imageSize will be done.
Convolver(const Convolver<FType>& other)
Convolver<FType> & operator=(const Convolver<FType> & other)
~Convolver()
void linearConv(Array<FType>& result, const Array<FType>& model, Bool fullSize=False)
void circularConv(Array<FType>& result, const Array<FType>& model)
void setPsf(const Array<FType>& psf, Bool cachePsf=False)
void setPsf(const Array<FType>& psf, IPosition imageShape, Bool fullSize=False, Bool cachePsf=False)
Set the transfer function for future convolutions to psf.
Assume linear convolution with a model of size imageSize
const Array<FType> getPsf(Bool cachePsf=True)
Get the psf currently used by this convolver
void setFastConvolve()
void makeXfr(const Array<FType>& psf, const IPosition& imageSize, Bool linear, Bool fullSize)
void makePsf(Array<FType>& psf)
IPosition defaultShape(const Array<FType>& psf)
IPosition extractShape(IPosition& psfSize, const IPosition& imageSize)
void doConvolution(Array<FType>& result, const Array<FType>& model, Bool fullSize)
void resizeXfr(const IPosition& imageShape, Bool linear, Bool fullSize)
void validate()