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

rfftw.hxx
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
60cosineTransformLineImpl(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
229realFourierTransformXEvenYEven(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
240realFourierTransformXEvenYEven(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::TinyVectorBase::begin
iterator begin()
Definition tinyvector.hxx:861
vigra
vigra::norm
FFTWComplex< R >::NormType norm(const FFTWComplex< R > &a)
norm (= magnitude)
Definition fftw3.hxx:1037

© 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.12.1