casacore
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 
39 namespace casacore { //# NAMESPACE CASACORE - BEGIN
40 
41 // Forward Declarations
42 template <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 
231 template<class FType> class Convolver
232 {
233 public:
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>
269  void linearConv(Array<FType>& result,
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 
303 private:
310 
311  void makeXfr(const Array<FType>& psf, const IPosition& imageSize,
312  Bool linear, Bool fullSize);
313  void makePsf(Array<FType>& psf);
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.
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
Convolver< FType > & operator=(const Convolver< FType > &other)
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:42
const Bool True
Definition: aipstype.h:43