[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

rfftw.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 1998-2002 by Ullrich Koethe */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_RFFTW_HXX
37 #define VIGRA_RFFTW_HXX
38 
39 #include "array_vector.hxx"
40 #include "fftw.hxx"
41 #include <rfftw.h>
42 
43 namespace vigra {
44 
45 namespace detail {
46 
47 struct FFTWSinCosConfig
48 {
49  ArrayVector<fftw_real> twiddles;
50  ArrayVector<fftw_real> fftwInput;
51  ArrayVector<fftw_complex> fftwTmpResult;
52  fftw_real norm;
53  rfftwnd_plan fftwPlan;
54 };
55 
56 template <class SrcIterator, class SrcAccessor,
57  class DestIterator, class DestAccessor,
58  class Config>
59 void
60 cosineTransformLineImpl(SrcIterator s, SrcIterator send, SrcAccessor src,
61  DestIterator d, DestAccessor dest,
62  Config & config)
63 {
64  int size = send - s;
65  int ns2 = size / 2;
66  int nm1 = size - 1;
67  int modn = size % 2;
68 
69  if(size <= 0)
70  return;
71 
72  fftw_real const * twiddles = config.twiddles.begin();
73  fftw_real * fftwInput = config.fftwInput.begin();
74  fftw_complex * fftwTmpResult = config.fftwTmpResult.begin();
75  fftw_real norm = config.norm;
76  rfftwnd_plan fftwPlan = config.fftwPlan;
77 
78  switch(size)
79  {
80  case 1:
81  {
82  dest.set(src(s) / norm, d);
83  break;
84  }
85  case 2:
86  {
87  dest.set((src(s) + src(s, 1)) / norm, d);
88  dest.set((src(s) - src(s, 1)) / norm, d, 1);
89  break;
90  }
91  case 3:
92  {
93  fftw_real x1p3 = src(s) + src(s, 2);
94  fftw_real tx2 = 2.0 * src(s, 1);
95 
96  dest.set((x1p3 + tx2) / norm, d, 0);
97  dest.set((src(s) - src(s, 2)) / norm, d, 1);
98  dest.set((x1p3 - tx2) / norm, d, 2);
99  break;
100  }
101  default:
102  {
103  fftw_real c1 = src(s) - src(s, nm1);
104  fftwInput[0] = src(s) + src(s, nm1);
105  for(int k=1; k<ns2; ++k)
106  {
107  int kc = nm1 - k;
108  fftw_real t1 = src(s, k) + src(s, kc);
109  fftw_real t2 = src(s, k) - src(s, kc);
110  c1 = c1 + twiddles[kc] * t2;
111  t2 = twiddles[k] * t2;
112  fftwInput[k] = t1 - t2;
113  fftwInput[kc] = t1 + t2;
114  }
115 
116  if (modn != 0)
117  {
118  fftwInput[ns2] = 2.0*src(s, ns2);
119  }
120  rfftwnd_one_real_to_complex(fftwPlan, fftwInput, fftwTmpResult);
121  dest.set(fftwTmpResult[0].re / norm, d, 0);
122  for(int k=1; k<(size+1)/2; ++k)
123  {
124  dest.set(fftwTmpResult[k].re, d, 2*k-1);
125  dest.set(fftwTmpResult[k].im, d, 2*k);
126  }
127  fftw_real xim2 = dest(d, 1);
128  dest.set(c1 / norm, d, 1);
129  for(int k=3; k<size; k+=2)
130  {
131  fftw_real xi = dest(d, k);
132  dest.set(dest(d, k-2) - dest(d, k-1) / norm, d, k);
133  dest.set(xim2 / norm, d, k-1);
134  xim2 = xi;
135  }
136  if (modn != 0)
137  {
138  dest.set(xim2 / norm, d, size-1);
139  }
140  }
141  }
142 }
143 
144 } // namespace detail
145 
146 template <class SrcTraverser, class SrcAccessor,
147  class DestTraverser, class DestAccessor>
148 void cosineTransformX(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
149  DestTraverser dul, DestAccessor dest, fftw_real norm)
150 {
151  int w = slr.x - sul.x;
152  int h = slr.y - sul.y;
153 
154  detail::FFTWSinCosConfig config;
155 
156  // horizontal transformation
157  int ns2 = w / 2;
158  int nm1 = w - 1;
159  int modn = w % 2;
160 
161  config.twiddles.resize(w+1);
162  config.fftwInput.resize(w+1);
163  config.fftwTmpResult.resize(w+1);
164  config.norm = norm;
165  config.fftwPlan = rfftw2d_create_plan(1, nm1, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE );
166 
167  fftw_real dt = M_PI / nm1;
168  for(int k=1; k<ns2; ++k)
169  {
170  fftw_real f = dt * k;
171  config.twiddles[k] = 2.0*VIGRA_CSTD::sin(f);
172  config.twiddles[nm1-k] = 2.0*VIGRA_CSTD::cos(f);
173  }
174 
175  for(; sul.y != slr.y; ++sul.y, ++dul.y)
176  {
177  typename SrcTraverser::row_iterator s = sul.rowIterator();
178  typename SrcTraverser::row_iterator send = s + w;
179  typename DestTraverser::row_iterator d = dul.rowIterator();
180  cosineTransformLineImpl(s, send, src, d, dest, config);
181  }
182 
183  rfftwnd_destroy_plan(config.fftwPlan);
184 }
185 
186 template <class SrcTraverser, class SrcAccessor,
187  class DestTraverser, class DestAccessor>
188 void cosineTransformY(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
189  DestTraverser dul, DestAccessor dest, fftw_real norm)
190 {
191  int w = slr.x - sul.x;
192  int h = slr.y - sul.y;
193 
194  detail::FFTWSinCosConfig config;
195 
196  // horizontal transformation
197  int ns2 = h / 2;
198  int nm1 = h - 1;
199  int modn = h % 2;
200 
201  config.twiddles.resize(h + 1);
202  config.fftwInput.resize(h + 1);
203  config.fftwTmpResult.resize(h + 1);
204  config.norm = norm;
205  config.fftwPlan = rfftw2d_create_plan(1, nm1, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE );
206 
207  fftw_real dt = M_PI / nm1;
208  for(int k=1; k<ns2; ++k)
209  {
210  fftw_real f = dt * k;
211  config.twiddles[k] = 2.0*VIGRA_CSTD::sin(f);
212  config.twiddles[nm1-k] = 2.0*VIGRA_CSTD::cos(f);
213  }
214 
215  for(; sul.x != slr.x; ++sul.x, ++dul.x)
216  {
217  typename SrcTraverser::column_iterator s = sul.columnIterator();
218  typename SrcTraverser::column_iterator send = s + h;
219  typename DestTraverser::column_iterator d = dul.columnIterator();
220  cosineTransformLineImpl(s, send, src, d, dest, config);
221  }
222 
223  rfftwnd_destroy_plan(config.fftwPlan);
224 }
225 
226 template <class SrcTraverser, class SrcAccessor,
227  class DestTraverser, class DestAccessor>
228 inline void
229 realFourierTransformXEvenYEven(SrcTraverser sul, SrcTraverser slr, SrcAccessor src,
230  DestTraverser dul, DestAccessor dest, fftw_real norm)
231 {
232  BasicImage<fftw_real> tmp(slr - sul);
233  cosineTransformX(sul, slr, src, tmp.upperLeft(), tmp.accessor(), 1.0);
234  cosineTransformY(tmp.upperLeft(), tmp.lowerRight(), tmp.accessor(), dul, dest, norm);
235 }
236 
237 template <class SrcTraverser, class SrcAccessor,
238  class DestTraverser, class DestAccessor>
239 inline void
240 realFourierTransformXEvenYEven(triple<SrcTraverser, SrcTraverser, SrcAccessor> src,
241  pair<DestTraverser, DestAccessor> dest, fftw_real norm)
242 {
243  realFourierTransformXEvenYEven(src.first, src.second, src.third, dest.first, dest.second, norm);
244 }
245 
246 } // namespace vigra
247 
248 #endif // VIGRA_RFFTW_HXX
vigra::linalg::sin
linalg::TemporaryMatrix< T > sin(MultiArrayView< 2, T, C > const &v)
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition: fftw3.hxx:1037
vigra::linalg::cos
linalg::TemporaryMatrix< T > cos(MultiArrayView< 2, T, C > const &v)

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)