npstat is hosted by Hepforge, IPPP Durham
NPStat  5.10.0
ConvolutionEngineND.hh
Go to the documentation of this file.
1 #ifndef NPSTAT_CONVOLUTIONENGINEND_HH_
2 #define NPSTAT_CONVOLUTIONENGINEND_HH_
3 
4 /*!
5 // \file ConvolutionEngineND.hh
6 //
7 // \brief Fast multidimensional convolutions via Fourier
8 // transforms (FFTW interface)
9 //
10 // Author: I. Volobouev
11 //
12 // December 2011
13 */
14 
15 #include <map>
16 #include <vector>
17 
19 
20 #ifdef SWIG
21 #include "npstat/wrap/arrayNDToNumpy.hh"
22 #include "npstat/wrap/numpyArrayUtils.hh"
23 #endif // SWIG
24 
25 namespace npstat {
26  /**
27  // Class for performing multivariate convolutions by FFT (using FFTW
28  // library). Internally, transforms are calculated in double precision.
29  */
31  {
32  public:
33  /**
34  // Constructor takes the shape of the data and FFTW optimization
35  // flag as parameters. Make sure the data dimensions are sufficient
36  // to avoid the convolution wrap-around.
37  */
38  ConvolutionEngineND(const unsigned *shape, unsigned rank,
39  unsigned optimization = FFTW_ESTIMATE);
41 
42  /** Provide the filter for subsequent convolutions */
43  template <typename Real>
44  void setFilter(const Real* data, const unsigned *dataShape,
45  unsigned dataRank);
46 
47  /**
48  // Convolve provided data with the filter previously established
49  // by setFilter
50  */
51  template <typename Real1, typename Real2>
52  void convolveWithFilter(const Real1* in, Real2* out,
53  const unsigned *dataShape, unsigned dataRank);
54 
55  /**
56  // Provide a filter for subsequent convolutions.
57  //
58  // "id" is the user id for the filter which can be later reused
59  // with the method "convolveWithDeposit".
60  */
61  template <typename Real>
62  void depositFilter(unsigned long id,
63  const Real* data, const unsigned *dataShape,
64  unsigned dataRank);
65 
66  /**
67  // Convolve provided data with one of the filters previously
68  // established by the depositFilter method. The "id" argument
69  // is used to identify a particular filter.
70  */
71  template <typename Real1, typename Real2>
72  void convolveWithDeposit(unsigned long id,
73  const Real1* in, Real2* out,
74  const unsigned *dataShape, unsigned dataRank);
75 
76  /**
77  // Convolve the filter established by the "setFilter" method
78  // with a previously deposited filter
79  */
80  template <typename Real2>
81  void convolveFilterAndDeposit(unsigned long id, Real2* out,
82  const unsigned *dataShape,
83  unsigned dataRank);
84 
85  /**
86  // Deconvolve the filter established by the "setFilter" method
87  // with a previously deposited filter
88  */
89  template <typename Real2>
90  void deconvolveFilterAndDeposit(unsigned long id, Real2* out,
91  const unsigned *dataShape,
92  unsigned dataRank);
93 
94  /**
95  // Discard a filter previously deposited by depositFilter method.
96  // This method returns "true" if the filter with the given id
97  // was indeed discarded and "false" if such a filter was not found.
98  */
99  bool discardFilter(unsigned long id);
100 
101  /** Discard all filters deposited by depositFilter method */
103 
104  /** Check if an array shape is compatible with this filter */
105  bool isShapeCompatible(const unsigned *shape, unsigned rank) const;
106 
107  /** Expected data dimensionality */
108  inline unsigned rank() const {return shape_.size();}
109 
110  /** Expected shape of the data */
111  inline const std::vector<unsigned>& shape() const {return shape_;}
112 
113  /** Expected data length */
114  inline unsigned long dataLength() const {return dataLen_;}
115 
116  // Some accessors
117  inline bool isFilterSet() const {return validFilter_;}
118  inline unsigned filterBankSize() const {return filterMap_.size();}
119 
120  private:
121  typedef std::map<unsigned long,FourierImage*> FilterMap;
122 
124  ConvolutionEngineND(const ConvolutionEngineND&);
125  ConvolutionEngineND& operator=(const ConvolutionEngineND&);
126 
127  template <typename Real1, typename Real2>
128  void convolveWithImage(const fftw_complex* image,
129  const Real1* in, Real2* out,
130  const unsigned *dataShape, unsigned dataRank);
131  fftw_plan pf;
132  fftw_plan pb;
133 
134  double* in;
135  fftw_complex* out;
136  fftw_complex* filterImage;
137 
138  std::vector<unsigned> shape_;
139  unsigned long dataLen_;
140  unsigned long cmplLen_;
141  bool validFilter_;
142 
143  FilterMap filterMap_;
144 
145 #ifdef SWIG
146  public:
147  inline void setFilter2(PyObject* input)
148  {
149  if (!isNumpyShapeCompatible(input, &shape_[0], shape_.size())) throw std::invalid_argument(
150  "In npstat::ConvolutionEngineND::setFilter2: incompatible array shape");
151  const int typenum = numpyArrayType(input);
152  if (typenum == NumpyTypecode<double>::code)
153  setFilter((double*)PyArray_DATA((PyArrayObject*)input), &shape_[0], shape_.size());
154  else if (typenum == NumpyTypecode<float>::code)
155  setFilter((float*)PyArray_DATA((PyArrayObject*)input), &shape_[0], shape_.size());
156  else if (typenum == NumpyTypecode<int>::code)
157  setFilter((int*)PyArray_DATA((PyArrayObject*)input), &shape_[0], shape_.size());
158  else if (typenum == NumpyTypecode<long long>::code)
159  setFilter((long long*)PyArray_DATA((PyArrayObject*)input), &shape_[0], shape_.size());
160  else if (typenum == NumpyTypecode<long int>::code)
161  setFilter((long int*)PyArray_DATA((PyArrayObject*)input), &shape_[0], shape_.size());
162  else if (typenum == NumpyTypecode<unsigned char>::code)
163  setFilter((unsigned char*)PyArray_DATA((PyArrayObject*)input), &shape_[0], shape_.size());
164  else
165  throw std::invalid_argument("In npstat::ConvolutionEngineND::setFilter2: "
166  "unsupported input array type");
167  }
168 
169  inline void depositFilter2(const unsigned long id,
170  PyObject* input)
171  {
172  const int typenum = numpyArrayType(input);
173  if (typenum != NumpyTypecode<double>::code) throw std::invalid_argument(
174  "In npstat::ConvolutionEngineND::depositFilter2: input is not an array of doubles");
175  if (!isNumpyShapeCompatible(input, &shape_[0], shape_.size())) throw std::invalid_argument(
176  "In npstat::ConvolutionEngineND::depositFilter2: incompatible array shape");
177  double* in = (double*)PyArray_DATA((PyArrayObject*)input);
178  depositFilter(id, in, &shape_[0], shape_.size());
179  }
180 
181  inline PyObject* convolveWithFilter2(PyObject* input,
182  const bool makeNonNegative = false)
183  {
184  const int typenum = numpyArrayType(input);
185  if (typenum != NumpyTypecode<double>::code) throw std::invalid_argument(
186  "In npstat::ConvolutionEngineND::convolveWithFilter2: input is not an array of doubles");
187  const unsigned dataRank = shape_.size();
188  assert(dataRank <= CHAR_BIT*sizeof(unsigned long));
189  if (!isNumpyShapeCompatible(input, &shape_[0], dataRank)) throw std::invalid_argument(
190  "In npstat::ConvolutionEngineND::convolveWithFilter2: incompatible array shape");
191  double* in = (double*)PyArray_DATA((PyArrayObject*)input);
192  npy_intp sh[CHAR_BIT*sizeof(unsigned long)];
193  for (unsigned i=0; i<dataRank; ++i)
194  sh[i] = shape_[i];
195  PyObject* xarr = PyArray_SimpleNew(dataRank, sh, typenum);
196  if (xarr)
197  {
198  double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
199  try
200  {
201  convolveWithFilter(in, out, &shape_[0], dataRank);
202  }
203  catch (const std::exception& e)
204  {
205  Py_DECREF(xarr);
206  throw e;
207  }
208  if (makeNonNegative)
209  for (unsigned long i=0; i<dataLen_; ++i)
210  if (out[i] < 0.0)
211  out[i] = 0.0;
212  }
213  return xarr;
214  }
215 
216  inline PyObject* convolveWithDeposit2(const unsigned long id,
217  PyObject* input,
218  const bool makeNonNegative = false)
219  {
220  const int typenum = numpyArrayType(input);
221  if (typenum != NumpyTypecode<double>::code) throw std::invalid_argument(
222  "In npstat::ConvolutionEngineND::convolveWithDeposit2: input is not an array of doubles");
223  const unsigned dataRank = shape_.size();
224  assert(dataRank <= CHAR_BIT*sizeof(unsigned long));
225  if (!isNumpyShapeCompatible(input, &shape_[0], dataRank)) throw std::invalid_argument(
226  "In npstat::ConvolutionEngineND::convolveWithDeposit2: incompatible array shape");
227  double* in = (double*)PyArray_DATA((PyArrayObject*)input);
228  npy_intp sh[CHAR_BIT*sizeof(unsigned long)];
229  for (unsigned i=0; i<dataRank; ++i)
230  sh[i] = shape_[i];
231  PyObject* xarr = PyArray_SimpleNew(dataRank, sh, typenum);
232  if (xarr)
233  {
234  double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
235  try
236  {
237  convolveWithDeposit(id, in, out, &shape_[0], dataRank);
238  }
239  catch (const std::exception& e)
240  {
241  Py_DECREF(xarr);
242  throw e;
243  }
244  if (makeNonNegative)
245  for (unsigned long i=0; i<dataLen_; ++i)
246  if (out[i] < 0.0)
247  out[i] = 0.0;
248  }
249  return xarr;
250  }
251 
252  inline PyObject* convolveFilterAndDeposit2(const unsigned long id,
253  const bool makeNonNegative = false)
254  {
255  const int typenum = NumpyTypecode<double>::code;
256  const unsigned dataRank = shape_.size();
257  assert(dataRank <= CHAR_BIT*sizeof(unsigned long));
258  npy_intp sh[CHAR_BIT*sizeof(unsigned long)];
259  for (unsigned i=0; i<dataRank; ++i)
260  sh[i] = shape_[i];
261  PyObject* xarr = PyArray_SimpleNew(dataRank, sh, typenum);
262  if (xarr)
263  {
264  double* out = (double*)PyArray_DATA((PyArrayObject*)xarr);
265  try
266  {
267  convolveFilterAndDeposit(id, out, &shape_[0], dataRank);
268  }
269  catch (const std::exception& e)
270  {
271  Py_DECREF(xarr);
272  throw e;
273  }
274  if (makeNonNegative)
275  for (unsigned long i=0; i<dataLen_; ++i)
276  if (out[i] < 0.0)
277  out[i] = 0.0;
278  }
279  return xarr;
280  }
281 #endif // SWIG
282  };
283 }
284 
285 #include "npstat/nm/ConvolutionEngineND.icc"
286 
287 #endif // NPSTAT_CONVOLUTIONENGINEND_HH_
Memory manager for Fourier transforms produced by FFTW.
Definition: ConvolutionEngineND.hh:31
bool discardFilter(unsigned long id)
void convolveFilterAndDeposit(unsigned long id, Real2 *out, const unsigned *dataShape, unsigned dataRank)
bool isShapeCompatible(const unsigned *shape, unsigned rank) const
unsigned long dataLength() const
Definition: ConvolutionEngineND.hh:114
void deconvolveFilterAndDeposit(unsigned long id, Real2 *out, const unsigned *dataShape, unsigned dataRank)
const std::vector< unsigned > & shape() const
Definition: ConvolutionEngineND.hh:111
void depositFilter(unsigned long id, const Real *data, const unsigned *dataShape, unsigned dataRank)
unsigned rank() const
Definition: ConvolutionEngineND.hh:108
void setFilter(const Real *data, const unsigned *dataShape, unsigned dataRank)
void convolveWithFilter(const Real1 *in, Real2 *out, const unsigned *dataShape, unsigned dataRank)
ConvolutionEngineND(const unsigned *shape, unsigned rank, unsigned optimization=FFTW_ESTIMATE)
void convolveWithDeposit(unsigned long id, const Real1 *in, Real2 *out, const unsigned *dataShape, unsigned dataRank)
Definition: AbsArrayProjector.hh:14