casacore
Loading...
Searching...
No Matches
Convolver.h
Go to the documentation of this file.
1//# Convolver.h: this defines Convolver a class for doing convolution
2//# Copyright (C) 1996,1999,2001,2002,2003
3//# Associated Universities, Inc. Washington DC, USA.
4//#
5//# This library is free software; you can redistribute it and/or modify it
6//# under the terms of the GNU Library General Public License as published by
7//# the Free Software Foundation; either version 2 of the License, or (at your
8//# option) any later version.
9//#
10//# This library is distributed in the hope that it will be useful, but WITHOUT
11//# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
12//# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
13//# License for more details.
14//#
15//# You should have received a copy of the GNU Library General Public License
16//# along with this library; if not, write to the Free Software Foundation,
17//# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
18//#
19//# Correspondence concerning AIPS++ should be addressed as follows:
20//# Internet email: aips2-request@nrao.edu.
21//# Postal address: AIPS++ Project Office
22//# National Radio Astronomy Observatory
23//# 520 Edgemont Road
24//# Charlottesville, VA 22903-2475 USA
25//#
26//#
27//# $Id$
28
29#ifndef SCIMATH_CONVOLVER_H
30#define SCIMATH_CONVOLVER_H
31
32#include <casacore/casa/aips.h>
33#include <casacore/casa/Arrays/Array.h>
34#include <casacore/casa/BasicSL/Complex.h>
35#include <casacore/scimath/Mathematics/FFTServer.h>
36#include <casacore/casa/Arrays/IPosition.h>
37#include <casacore/scimath/Mathematics/NumericTraits.h>
38
39namespace casacore { //# NAMESPACE CASACORE - BEGIN
40
41// Forward Declarations
42template <class FType> class Convolver;
43
44// Typedefs
47
48// <summary>
49// A class for doing multi-dimensional convolution
50// </summary>
51
52// <use visibility=export>
53
54// <reviewed reviewer="Brian Glendenning " date="1996/05/27"
55// tests="tConvolution">
56// </reviewed>
57
58// <prerequisite>
59// <li> The mathematical concept of convolution
60// </prerequisite>
61//
62// <etymology>
63// The convolver class performs convolution!
64// </etymology>
65//
66// <synopsis>
67// This class will perform linear or circular convolution on arrays.
68//
69// The dimension of the convolution done is determined by the dimension of
70// the point spread function (psf), so for example, if the psf is a Vector,
71// one dimensional convolution will be done. The dimension of the model
72// that is to be convolved must be at least the same as the point
73// spread function, but it can be larger. If it is then the convolution will
74// be repeated for each row or plane of the model.
75// <note role=tip>
76// This class strips all degenerate axes when determining the dimensionality
77// of the psf or model. So a psf with shapes of [1,1,16] or [16,1,1] is
78// treated the same as a Vector of length 16, and will result in one
79// dimensional convolutions along the first non-degenerate axis of the
80// supplied model.
81// </note>
82
83// Repeated convolution can only be done along the fastest moving axes of
84// the supplied image. For example, if a one dimensional psf is used (so
85// that one dimensional convolution is being done), and a cube of data is
86// supplied then the convolution will be repeated for each row in the
87// cube. It is not currently possible to have this class do repeated one
88// dimensional convolution along all the columns or along the z
89// axis. To do this you need to use an iterator external to the class to
90// successively feed in the appropriate slices of your Array.
91
92// The difference between linear and circular convolution can best be
93// explained with a one dimensional example.
94// Suppose the psf and data to be convolved are:
95// <srcblock>
96// psf = [0 .5 1 .1]; data = [1 0 0 0 0 0]
97// </srcblock>
98// then their linear and circular convolutions are:
99// <srcblock>
100// circular convolution = [1 .1 0 0 0 .5]
101// linear convolution = [1 .1 0 0 0 0] (fullSize == False)
102// linear convolution = [0 .5 1 .1 0 0 0 0 0] (fullSize == True)
103// </srcblock>
104// The circular convolution "wraps around" whereas the linear one does not.
105// Usage of the fullSize option is explained below. As can be seen from the
106// above example this class does not normalise the convolved result by any
107// factor that depends on the psf, so if the "beam area" is not unity the
108// flux scales will vary.
109
110// The "centre" of the convolution is at the point (NX/2, NY/2) (assuming a
111// 2 dimensional psf) where the first point in the psf is at (0,0) and the
112// last is at (NX-1, NY-1). This means that a psf that is all zero except
113// for 1 at the "centre" pixel will when convolved with any model leave the
114// model unchanged.
115
116// The convolution is done in the Fourier domain and the transform of the
117// psf (the transfer function) is cached by this class. If the cached
118// transfer function is the wrong size for a given model it will be
119// automatically be recomputed to the right size (this will involve two
120// FFT's)
121
122// Each convolution requires two Fourier transforms which dominate the
123// computational load. Hence the computational expense is
124// <em> n Log(n) </em> for 1 dimensional and
125// <em> n^2 Log(n) </em> for 2 dimensional convolutions.
126
127// The size of the convolved result is always the same as the input model
128// unless linear convolution is done with the fullSize option set to True.
129// In this case the result will be larger than the model and include the
130// full linear convolution (resultSize = psfSize+modelSize-1), rather than
131// the central portion.
132
133// If the convolver is constructed with an expected model size (as in the
134// example below) then the cached transfer function will be computed to a
135// size appropriate for linear convolution of models of that size. If no
136// model size is given then the cached transfer function will be computed
137// with a size appropriate for circular convolution. These guidelines also
138// apply when using the setPsf functions.
139
140// <note role=tip>
141// If you are intending to do 'fullsize' linear convolutions
142// you should also set the fullsize option to True as the cached transfer
143// function is a different size for fullsize linear convolutions.
144// </note>
145
146// For linear convolution the psf can be larger, the same size or smaller
147// than the model but for circular convolution the psf must be smaller or the
148// same size.
149
150// The size of the cached transfer function (and also the length of the
151// FFT's calculated) depends on the sizes of the psf and the model, as well
152// as whether you are doing linear or circular convolution and the fullSize
153// option. It is always advantageous to use the smallest possible psf
154// (ie. do not pad the psf prior to supplying it to this class). Be aware
155// that using odd length images will lead to this class doing odd length
156// FFT's, which are less computationally effecient (particularly is the
157// length of the transform is a prime number) in general than even length
158// transforms.
159
160// There are only two valid template types namely,
161// <ol>
162// <li>FType=Float or
163// <li>FType=Double
164// </ol>
165// and the user may prefer to use the following typedef's:
166// <srcblock>
167// FloatConvolver (= Convolver<Float>) or
168// DoubleConvolver (= Convolver<Double>)
169// </srcblock>
170// rather than explicitly specifying the template arguements.
171// <note role=tip>
172// The typedefs need to be redeclared when using the gnu compiler making
173// them essentially useless.
174// </note>
175
176// When this class is constructed you may choose to have the psf
177// explicitly stored by the class (by setting cachePsf=True). This will
178// allow speedy access to the psf when using the getPsf function. However
179// the getPsf function can still be called even if the psf has not been
180// cached. Then the psf will be computed by FFT'ing the transfer
181// function, and the psf will also then be cached (unless
182// cachePsf=Flase). Cacheing the psf is also a good idea if you will be
183// switching between different sized transfer functions (eg. mixing
184// linear and circular convolution) as it will save one of the two
185// FFT's. Note that even though the psf is returned as a const Array, it
186// is possible to inadvertently modify it using the array copy constructor
187// as this uses reference symantics. Modifying the psf is NOT
188// recommended. eg.
189// <srcblock>
190// DoubleConvolver conv();
191// {
192// Matrix<Double> psf(20,20);
193// conv.setPsf(psf);
194// }
195// Matrix<Double> convPsf = conv.getPsf(); // Get the psf used by the convolver
196// convPsf(0,0) = -100; // And modify it. This modifies
197// // This internal psf used by the
198// // convolver also! (unless it is
199// // not caching the psf)
200// </srcblock>
201//
202// </synopsis>
203//
204// <example>
205// Calculate the convolution of two Matrices (psf and model);
206// <srcblock>
207// Matrix<Float> psf(4,4), model(12,12);
208// ...put meaningful values into the above two matrices...
209// FloatConvolver conv(psf, model.shape());
210// conv.linearConv(result, model); // result = Convolution(psf, model)
211// </srcblock>
212// </example>
213//
214// <motivation>
215// I needed to do linear convolution to write a clean algorithm. It
216// blossomed into this class.
217// </motivation>
218//
219// <thrown>
220// <li> AipsError: if psf has more dimensions than the model.
221// </thrown>
222//
223// <todo asof="yyyy/mm/dd">
224// <li> the class should detect if the psf or image is small and do the
225// convolution directly rather than use the Fourier domain
226// <li> Add a lattice interface, and more flexible iteration scheme
227// <li> Allow the psf to be specified with a
228// <linkto class=Function>Function</linkto>.
229// </todo>
230
231template<class FType> class Convolver
232{
233public:
234 // When using the default constructor the psf MUST be specified using the
235 // setPsf function prior to doing any convolution.
236 // <group>
238 // </group>
239 // Create the cached Transfer function assuming that circular convolution
240 // will be done
241 // <group>
242 Convolver(const Array<FType>& psf, Bool cachePsf=False);
243 // </group>
244 // Create the cached Transfer function assuming that linear convolution
245 // with an array of size imageSize will be done.
246 // <group>
247 Convolver(const Array<FType>& psf, const IPosition& imageSize,
248 Bool fullSize=False, Bool cachePsf=False);
249 // </group>
250
251 // The copy constructor and the assignment operator make copies (and not
252 // references) of all the internal data arrays, as this object could get
253 // really screwed up if the private data was silently messed with.
254 // <group>
257 // </group>
258
259 // The destructor does nothing!
260 // <group>
262 // </group>
263
264 // Perform linear convolution of the model with the previously
265 // specified psf. Return the answer in result. Set fullSize to True if you
266 // want the full convolution, rather than the central portion (the same
267 // size as the model) returned.
268 // <group>
270 const Array<FType>& model,
271 Bool fullSize=False);
272 // </group>
273
274 // Perform circular convolution of the model with the previously
275 // specified psf. Return the answer in result.
276 // <group>
278 const Array<FType>& model);
279 // </group>
280
281 // Set the transfer function for future convolutions to psf.
282 // Assume circular convolution will be done
283 // <group>
284 void setPsf(const Array<FType>& psf, Bool cachePsf=False);
285 // </group>
286 // Set the transfer function for future convolutions to psf.
287 // Assume linear convolution with a model of size imageSize
288 // <group>
289 void setPsf(const Array<FType>& psf,
290 IPosition imageShape, Bool fullSize=False,
291 Bool cachePsf=False);
292 // </group>
293 // Get the psf currently used by this convolver
294 // <group>
295 const Array<FType> getPsf(Bool cachePsf=True);
296 // </group>
297
298 // Set to use convolution with lesser flips
299 // <group>
301 // </group>
302
303private:
310
311 void makeXfr(const Array<FType>& psf, const IPosition& imageSize,
312 Bool linear, Bool fullSize);
315 IPosition extractShape(IPosition& psfSize, const IPosition& imageSize);
317 const Array<FType>& model,
318 Bool fullSize);
319 void resizeXfr(const IPosition& imageShape, Bool linear, Bool fullSize);
320//# void padArray(Array<FType>& paddedArr, const Array<FType>& origArr,
321//# const IPosition & blc);
324 void validate();
325};
326
327} //# NAMESPACE CASACORE - END
328
329#ifndef CASACORE_NO_AUTO_TEMPLATES
330#include <casacore/scimath/Mathematics/Convolver.tcc>
331#endif //# CASACORE_NO_AUTO_TEMPLATES
332#endif
Forward Declarations.
Definition Convolver.h:232
Array< FType > thePsf
Definition Convolver.h:307
void setFastConvolve()
Set to use convolution with lesser flips.
Convolver< FType > & operator=(const Convolver< FType > &other)
void setPsf(const Array< FType > &psf, Bool cachePsf=False)
Set the transfer function for future convolutions to psf.
void makePsf(Array< FType > &psf)
FFTServer< FType, typename NumericTraits< FType >::ConjugateType > theIFFT
Definition Convolver.h:309
Convolver(const Array< FType > &psf, Bool cachePsf=False)
Create the cached Transfer function assuming that circular convolution will be done.
const Array< FType > getPsf(Bool cachePsf=True)
Get the psf currently used by this convolver.
void linearConv(Array< FType > &result, const Array< FType > &model, Bool fullSize=False)
Perform linear convolution of the model with the previously specified psf.
void circularConv(Array< FType > &result, const Array< FType > &model)
Perform circular convolution of the model with the previously specified psf.
void makeXfr(const Array< FType > &psf, const IPosition &imageSize, Bool linear, Bool fullSize)
Convolver(const Convolver< FType > &other)
The copy constructor and the assignment operator make copies (and not references) of all the internal...
void resizeXfr(const IPosition &imageShape, Bool linear, Bool fullSize)
IPosition extractShape(IPosition &psfSize, const IPosition &imageSize)
IPosition thePsfSize
Definition Convolver.h:304
~Convolver()
The destructor does nothing!
IPosition defaultShape(const Array< FType > &psf)
Convolver()
When using the default constructor the psf MUST be specified using the setPsf function prior to doing...
Definition Convolver.h:237
void setPsf(const Array< FType > &psf, IPosition imageShape, Bool fullSize=False, Bool cachePsf=False)
Set the transfer function for future convolutions to psf.
FFTServer< FType, typename NumericTraits< FType >::ConjugateType > theFFT
Definition Convolver.h:308
Array< typename NumericTraits< FType >::ConjugateType > theXfr
Definition Convolver.h:306
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 ...
void doConvolution(Array< FType > &result, const Array< FType > &model, Bool fullSize)
IPosition theFFTSize
Definition Convolver.h:305
A class with methods for Fast Fourier Transforms.
Definition FFTServer.h:228
this file contains all the compiler specific defines
Definition mainpage.dox:28
const Bool False
Definition aipstype.h:44
Convolver< Double > DoubleConvolver
Definition Convolver.h:46
bool Bool
Define the standard types used by Casacore.
Definition aipstype.h:42
Convolver< Float > FloatConvolver
Typedefs.
Definition Convolver.h:45
const Bool True
Definition aipstype.h:43