[ VIGRA Homepage | Class Index | Function Index | File Index | Main Page ]
![]() |
vigra/rfftw.hxx | ![]() |
---|
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2002 by Ullrich Koethe */ 00004 /* Cognitive Systems Group, University of Hamburg, Germany */ 00005 /* */ 00006 /* This file is part of the VIGRA computer vision library. */ 00007 /* ( Version 1.3.2, Jan 27 2005 ) */ 00008 /* You may use, modify, and distribute this software according */ 00009 /* to the terms stated in the LICENSE file included in */ 00010 /* the VIGRA distribution. */ 00011 /* */ 00012 /* The VIGRA Website is */ 00013 /* http://kogs-www.informatik.uni-hamburg.de/~koethe/vigra/ */ 00014 /* Please direct questions, bug reports, and contributions to */ 00015 /* koethe@informatik.uni-hamburg.de */ 00016 /* */ 00017 /* THIS SOFTWARE IS PROVIDED AS IS AND WITHOUT ANY EXPRESS OR */ 00018 /* IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED */ 00019 /* WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE. */ 00020 /* */ 00021 /************************************************************************/ 00022 00023 #ifndef VIGRA_RFFTW_HXX 00024 #define VIGRA_RFFTW_HXX 00025 00026 #include "vigra/array_vector.hxx" 00027 #include "vigra/fftw.hxx" 00028 #include <rfftw.h> 00029 00030 namespace vigra { 00031 00032 namespace detail { 00033 00034 struct FFTWSinCosConfig 00035 { 00036 ArrayVector<fftw_real> twiddles; 00037 ArrayVector<fftw_real> fftwInput; 00038 ArrayVector<fftw_complex> fftwTmpResult; 00039 fftw_real norm; 00040 rfftwnd_plan fftwPlan; 00041 }; 00042 00043 template <class SrcIterator, class SrcAccessor, 00044 class DestIterator, class DestAccessor, 00045 class Config> 00046 void 00047 cosineTransformLineImpl(SrcIterator s, SrcIterator send, SrcAccessor src, 00048 DestIterator d, DestAccessor dest, 00049 Config & config) 00050 { 00051 int size = send - s; 00052 int ns2 = size / 2; 00053 int nm1 = size - 1; 00054 int modn = size % 2; 00055 00056 if(size <= 0) 00057 return; 00058 00059 fftw_real const * twiddles = config.twiddles.begin(); 00060 fftw_real * fftwInput = config.fftwInput.begin(); 00061 fftw_complex * fftwTmpResult = config.fftwTmpResult.begin(); 00062 fftw_real norm = config.norm; 00063 rfftwnd_plan fftwPlan = config.fftwPlan; 00064 00065 switch(size) 00066 { 00067 case 1: 00068 { 00069 dest.set(src(s) / norm, d); 00070 break; 00071 } 00072 case 2: 00073 { 00074 dest.set((src(s) + src(s, 1)) / norm, d); 00075 dest.set((src(s) - src(s, 1)) / norm, d, 1); 00076 break; 00077 } 00078 case 3: 00079 { 00080 fftw_real x1p3 = src(s) + src(s, 2); 00081 fftw_real tx2 = 2.0 * src(s, 1); 00082 00083 dest.set((x1p3 + tx2) / norm, d, 0); 00084 dest.set((src(s) - src(s, 2)) / norm, d, 1); 00085 dest.set((x1p3 - tx2) / norm, d, 2); 00086 break; 00087 } 00088 default: 00089 { 00090 fftw_real c1 = src(s) - src(s, nm1); 00091 fftwInput[0] = src(s) + src(s, nm1); 00092 for(int k=1; k<ns2; ++k) 00093 { 00094 int kc = nm1 - k; 00095 fftw_real t1 = src(s, k) + src(s, kc); 00096 fftw_real t2 = src(s, k) - src(s, kc); 00097 c1 = c1 + twiddles[kc] * t2; 00098 t2 = twiddles[k] * t2; 00099 fftwInput[k] = t1 - t2; 00100 fftwInput[kc] = t1 + t2; 00101 } 00102 00103 if (modn != 0) 00104 { 00105 fftwInput[ns2] = 2.0*src(s, ns2); 00106 } 00107 rfftwnd_one_real_to_complex(fftwPlan, fftwInput, fftwTmpResult); 00108 dest.set(fftwTmpResult[0].re / norm, d, 0); 00109 for(int k=1; k<(size+1)/2; ++k) 00110 { 00111 dest.set(fftwTmpResult[k].re, d, 2*k-1); 00112 dest.set(fftwTmpResult[k].im, d, 2*k); 00113 } 00114 fftw_real xim2 = dest(d, 1); 00115 dest.set(c1 / norm, d, 1); 00116 for(int k=3; k<size; k+=2) 00117 { 00118 fftw_real xi = dest(d, k); 00119 dest.set(dest(d, k-2) - dest(d, k-1) / norm, d, k); 00120 dest.set(xim2 / norm, d, k-1); 00121 xim2 = xi; 00122 } 00123 if (modn != 0) 00124 { 00125 dest.set(xim2 / norm, d, size-1); 00126 } 00127 } 00128 } 00129 } 00130 00131 } // namespace detail 00132 00133 template <class SrcTraverser, class SrcAccessor, 00134 class DestTraverser, class DestAccessor> 00135 void cosineTransformX(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 00136 DestTraverser dul, DestAccessor dest, fftw_real norm) 00137 { 00138 int w = slr.x - sul.x; 00139 int h = slr.y - sul.y; 00140 00141 detail::FFTWSinCosConfig config; 00142 00143 // horizontal transformation 00144 int ns2 = w / 2; 00145 int nm1 = w - 1; 00146 int modn = w % 2; 00147 00148 config.twiddles.resize(w+1); 00149 config.fftwInput.resize(w+1); 00150 config.fftwTmpResult.resize(w+1); 00151 config.norm = norm; 00152 config.fftwPlan = rfftw2d_create_plan(1, nm1, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE ); 00153 00154 fftw_real dt = M_PI / nm1; 00155 for(int k=1; k<ns2; ++k) 00156 { 00157 fftw_real f = dt * k; 00158 config.twiddles[k] = 2.0*VIGRA_CSTD::sin(f); 00159 config.twiddles[nm1-k] = 2.0*VIGRA_CSTD::cos(f); 00160 } 00161 00162 for(; sul.y != slr.y; ++sul.y, ++dul.y) 00163 { 00164 typename SrcTraverser::row_iterator s = sul.rowIterator(); 00165 typename SrcTraverser::row_iterator send = s + w; 00166 typename DestTraverser::row_iterator d = dul.rowIterator(); 00167 cosineTransformLineImpl(s, send, src, d, dest, config); 00168 } 00169 00170 rfftwnd_destroy_plan(config.fftwPlan); 00171 } 00172 00173 template <class SrcTraverser, class SrcAccessor, 00174 class DestTraverser, class DestAccessor> 00175 void cosineTransformY(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 00176 DestTraverser dul, DestAccessor dest, fftw_real norm) 00177 { 00178 int w = slr.x - sul.x; 00179 int h = slr.y - sul.y; 00180 00181 detail::FFTWSinCosConfig config; 00182 00183 // horizontal transformation 00184 int ns2 = h / 2; 00185 int nm1 = h - 1; 00186 int modn = h % 2; 00187 00188 config.twiddles.resize(h + 1); 00189 config.fftwInput.resize(h + 1); 00190 config.fftwTmpResult.resize(h + 1); 00191 config.norm = norm; 00192 config.fftwPlan = rfftw2d_create_plan(1, nm1, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE ); 00193 00194 fftw_real dt = M_PI / nm1; 00195 for(int k=1; k<ns2; ++k) 00196 { 00197 fftw_real f = dt * k; 00198 config.twiddles[k] = 2.0*VIGRA_CSTD::sin(f); 00199 config.twiddles[nm1-k] = 2.0*VIGRA_CSTD::cos(f); 00200 } 00201 00202 for(; sul.x != slr.x; ++sul.x, ++dul.x) 00203 { 00204 typename SrcTraverser::column_iterator s = sul.columnIterator(); 00205 typename SrcTraverser::column_iterator send = s + h; 00206 typename DestTraverser::column_iterator d = dul.columnIterator(); 00207 cosineTransformLineImpl(s, send, src, d, dest, config); 00208 } 00209 00210 rfftwnd_destroy_plan(config.fftwPlan); 00211 } 00212 00213 template <class SrcTraverser, class SrcAccessor, 00214 class DestTraverser, class DestAccessor> 00215 inline void 00216 realFourierTransformXEvenYEven(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 00217 DestTraverser dul, DestAccessor dest, fftw_real norm) 00218 { 00219 BasicImage<fftw_real> tmp(slr - sul); 00220 cosineTransformX(sul, slr, src, tmp.upperLeft(), tmp.accessor(), 1.0); 00221 cosineTransformY(tmp.upperLeft(), tmp.lowerRight(), tmp.accessor(), dul, dest, norm); 00222 } 00223 00224 template <class SrcTraverser, class SrcAccessor, 00225 class DestTraverser, class DestAccessor> 00226 inline void 00227 realFourierTransformXEvenYEven(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 00228 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 00229 { 00230 realFourierTransformXEvenYEven(src.first, src.second, src.third, dest.first, dest.second, norm); 00231 } 00232 00233 } // namespace vigra 00234 00235 #endif // VIGRA_RFFTW_HXX
© Ullrich Köthe (koethe@informatik.uni-hamburg.de) |
html generated using doxygen and Python
|