Loading...
Searching...
No Matches
Convolver.h
Go to the documentation of this file.
//# Convolver.h: this defines Convolver a class for doing convolution
//# Copyright (C) 1996,1999,2001,2002,2003
//# Associated Universities, Inc. Washington DC, USA.
//#
//# This library is free software; you can redistribute it and/or modify it
//# under the terms of the GNU Library General Public License as published by
//# the Free Software Foundation; either version 2 of the License, or (at your
//# option) any later version.
//#
//# This library is distributed in the hope that it will be useful, but WITHOUT
//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
//# License for more details.
//#
//# You should have received a copy of the GNU Library General Public License
//# along with this library; if not, write to the Free Software Foundation,
//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
//#
//# Correspondence concerning AIPS++ should be addressed as follows:
//# Internet email: casa-feedback@nrao.edu.
//# Postal address: AIPS++ Project Office
//# National Radio Astronomy Observatory
//# 520 Edgemont Road
//# Charlottesville, VA 22903-2475 USA
#ifndef SCIMATH_CONVOLVER_H
#define SCIMATH_CONVOLVER_H
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Array.h>
#include <casacore/casa/BasicSL/Complex.h>
#include <casacore/scimath/Mathematics/FFTServer.h>
#include <casacore/casa/Arrays/IPosition.h>
#include <casacore/scimath/Mathematics/NumericTraits.h>
namespace casacore { //# NAMESPACE CASACORE - BEGIN
// Forward Declarations
template <class FType> class Convolver;
// Typedefs
typedef Convolver<Float> FloatConvolver;
typedef Convolver<Double> DoubleConvolver;
// <summary>
// A class for doing multi-dimensional convolution
// </summary>
// <use visibility=export>
// <reviewed reviewer="Brian Glendenning " date="1996/05/27"
// tests="tConvolution">
// </reviewed>
// <prerequisite>
// <li> The mathematical concept of convolution
// </prerequisite>
//
// <etymology>
// The convolver class performs convolution!
// </etymology>
//
// <synopsis>
// This class will perform linear or circular convolution on arrays.
//
// 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.
// <note role=tip>
// 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.
// </note>
// 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:
// <srcblock>
// psf = [0 .5 1 .1]; data = [1 0 0 0 0 0]
// </srcblock>
// then their linear and circular convolutions are:
// <srcblock>
// 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)
// </srcblock>
// 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.
// 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
// <em> n Log(n) </em> for 1 dimensional and
// <em> n^2 Log(n) </em> 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.
// <note role=tip>
// 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.
// </note>
// 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,
// <ol>
// <li>FType=Float or
// <li>FType=Double
// </ol>
// and the user may prefer to use the following typedef's:
// <srcblock>
// FloatConvolver (= Convolver<Float>) or
// DoubleConvolver (= Convolver<Double>)
// </srcblock>
// rather than explicitly specifying the template arguements.
// <note role=tip>
// The typedefs need to be redeclared when using the gnu compiler making
// them essentially useless.
// </note>
// 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.
// <srcblock>
// 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)
// </srcblock>
//
// </synopsis>
//
// <example>
// Calculate the convolution of two Matrices (psf and model);
// <srcblock>
// 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)
// </srcblock>
// </example>
//
// <motivation>
// I needed to do linear convolution to write a clean algorithm. It
// blossomed into this class.
// </motivation>
//
// <thrown>
// <li> AipsError: if psf has more dimensions than the model.
// </thrown>
//
// <todo asof="yyyy/mm/dd">
// <li> the class should detect if the psf or image is small and do the
// convolution directly rather than use the Fourier domain
// <li> Add a lattice interface, and more flexible iteration scheme
// <li> Allow the psf to be specified with a
// <linkto class=Function>Function</linkto>.
// </todo>
template<class FType> class Convolver
{
public:
// When using the default constructor the psf MUST be specified using the
// setPsf function prior to doing any convolution.
// <group>
Convolver(){}
// </group>
// Create the cached Transfer function assuming that circular convolution
// will be done
// <group>
Convolver(const Array<FType>& psf, Bool cachePsf=False);
// </group>
// Create the cached Transfer function assuming that linear convolution
// with an array of size imageSize will be done.
// <group>
Convolver(const Array<FType>& psf, const IPosition& imageSize,
Bool fullSize=False, Bool cachePsf=False);
// </group>
// 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.
// <group>
Convolver(const Convolver<FType>& other);
Convolver<FType> & operator=(const Convolver<FType> & other);
// </group>
// The destructor does nothing!
// <group>
~Convolver();
// </group>
// 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.
// <group>
void linearConv(Array<FType>& result,
const Array<FType>& model,
Bool fullSize=False);
// </group>
// Perform circular convolution of the model with the previously
// specified psf. Return the answer in result.
// <group>
void circularConv(Array<FType>& result,
const Array<FType>& model);
// </group>
// Set the transfer function for future convolutions to psf.
// Assume circular convolution will be done
// <group>
void setPsf(const Array<FType>& psf, Bool cachePsf=False);
// </group>
// Set the transfer function for future convolutions to psf.
// Assume linear convolution with a model of size imageSize
// <group>
void setPsf(const Array<FType>& psf,
IPosition imageShape, Bool fullSize=False,
Bool cachePsf=False);
// </group>
// Get the psf currently used by this convolver
// <group>
const Array<FType> getPsf(Bool cachePsf=True);
// </group>
// Set to use convolution with lesser flips
// <group>
void setFastConvolve();
// </group>
private:
IPosition thePsfSize;
IPosition theFFTSize;
Array<typename NumericTraits<FType>::ConjugateType> theXfr;
Array<FType> thePsf;
FFTServer<FType, typename NumericTraits<FType>::ConjugateType> theFFT;
FFTServer<FType, typename NumericTraits<FType>::ConjugateType> theIFFT;
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 padArray(Array<FType>& paddedArr, const Array<FType>& origArr,
//# const IPosition & blc);
Bool valid;
Bool doFast_p;
void validate();
};
} //# NAMESPACE CASACORE - END
#ifndef CASACORE_NO_AUTO_TEMPLATES
#include <casacore/scimath/Mathematics/Convolver.tcc>
#endif //# CASACORE_NO_AUTO_TEMPLATES
#endif
casacore::Array
Definition Array.h:155
casacore::Convolver
Forward Declarations.
Definition Convolver.h:229
casacore::Convolver< Float >::thePsf
Array< Float > thePsf
Definition Convolver.h:304
casacore::Convolver::setFastConvolve
void setFastConvolve()
Set to use convolution with lesser flips.
casacore::Convolver::operator=
Convolver< FType > & operator=(const Convolver< FType > &other)
casacore::Convolver::setPsf
void setPsf(const Array< FType > &psf, Bool cachePsf=False)
Set the transfer function for future convolutions to psf.
casacore::Convolver::makePsf
void makePsf(Array< FType > &psf)
casacore::Convolver< Float >::theIFFT
FFTServer< Float, typename NumericTraits< Float >::ConjugateType > theIFFT
Definition Convolver.h:306
casacore::Convolver::Convolver
Convolver(const Array< FType > &psf, Bool cachePsf=False)
Create the cached Transfer function assuming that circular convolution will be done.
casacore::Convolver::getPsf
const Array< FType > getPsf(Bool cachePsf=True)
Get the psf currently used by this convolver.
casacore::Convolver< Float >::valid
Bool valid
Definition Convolver.h:319
casacore::Convolver::linearConv
void linearConv(Array< FType > &result, const Array< FType > &model, Bool fullSize=False)
Perform linear convolution of the model with the previously specified psf.
casacore::Convolver::circularConv
void circularConv(Array< FType > &result, const Array< FType > &model)
Perform circular convolution of the model with the previously specified psf.
casacore::Convolver::makeXfr
void makeXfr(const Array< FType > &psf, const IPosition &imageSize, Bool linear, Bool fullSize)
casacore::Convolver::Convolver
Convolver(const Convolver< FType > &other)
The copy constructor and the assignment operator make copies (and not references) of all the internal...
casacore::Convolver::validate
void validate()
casacore::Convolver::resizeXfr
void resizeXfr(const IPosition &imageShape, Bool linear, Bool fullSize)
casacore::Convolver< Float >::doFast_p
Bool doFast_p
Definition Convolver.h:320
casacore::Convolver::extractShape
IPosition extractShape(IPosition &psfSize, const IPosition &imageSize)
casacore::Convolver< Float >::thePsfSize
IPosition thePsfSize
Definition Convolver.h:301
casacore::Convolver::~Convolver
~Convolver()
The destructor does nothing!
casacore::Convolver::defaultShape
IPosition defaultShape(const Array< FType > &psf)
casacore::Convolver::Convolver
Convolver()
When using the default constructor the psf MUST be specified using the setPsf function prior to doing...
Definition Convolver.h:234
casacore::Convolver::setPsf
void setPsf(const Array< FType > &psf, IPosition imageShape, Bool fullSize=False, Bool cachePsf=False)
Set the transfer function for future convolutions to psf.
casacore::Convolver< Float >::theFFT
FFTServer< Float, typename NumericTraits< Float >::ConjugateType > theFFT
Definition Convolver.h:305
casacore::Convolver< Float >::theXfr
Array< typename NumericTraits< Float >::ConjugateType > theXfr
Definition Convolver.h:303
casacore::Convolver::Convolver
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 ...
casacore::Convolver::doConvolution
void doConvolution(Array< FType > &result, const Array< FType > &model, Bool fullSize)
casacore::Convolver< Float >::theFFTSize
IPosition theFFTSize
Definition Convolver.h:302
casacore::FFTServer
A class with methods for Fast Fourier Transforms.
Definition FFTServer.h:226
casacore::IPosition
Definition IPosition.h:117
casacore
this file contains all the compiler specific defines
Definition mainpage.dox:28
casacore::False
const Bool False
Definition aipstype.h:42
casacore::DoubleConvolver
Convolver< Double > DoubleConvolver
Definition Convolver.h:43
casacore::Bool
bool Bool
Define the standard types used by Casacore.
Definition aipstype.h:40
casacore::FloatConvolver
Convolver< Float > FloatConvolver
Typedefs.
Definition Convolver.h:42
casacore::True
const Bool True
Definition aipstype.h:41