/usr/include/casacore/scimath/Mathematics/Convolver.h is in casacore-dev 2.2.0-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 | //# 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: aips2-request@nrao.edu.
//# Postal address: AIPS++ Project Office
//# National Radio Astronomy Observatory
//# 520 Edgemont Road
//# Charlottesville, VA 22903-2475 USA
//#
//#
//# $Id$
#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
|