[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]
vigra/rfftw.hxx | ![]() |
00001 /************************************************************************/ 00002 /* */ 00003 /* Copyright 1998-2002 by Ullrich Koethe */ 00004 /* */ 00005 /* This file is part of the VIGRA computer vision library. */ 00006 /* The VIGRA Website is */ 00007 /* http://hci.iwr.uni-heidelberg.de/vigra/ */ 00008 /* Please direct questions, bug reports, and contributions to */ 00009 /* ullrich.koethe@iwr.uni-heidelberg.de or */ 00010 /* vigra@informatik.uni-hamburg.de */ 00011 /* */ 00012 /* Permission is hereby granted, free of charge, to any person */ 00013 /* obtaining a copy of this software and associated documentation */ 00014 /* files (the "Software"), to deal in the Software without */ 00015 /* restriction, including without limitation the rights to use, */ 00016 /* copy, modify, merge, publish, distribute, sublicense, and/or */ 00017 /* sell copies of the Software, and to permit persons to whom the */ 00018 /* Software is furnished to do so, subject to the following */ 00019 /* conditions: */ 00020 /* */ 00021 /* The above copyright notice and this permission notice shall be */ 00022 /* included in all copies or substantial portions of the */ 00023 /* Software. */ 00024 /* */ 00025 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */ 00026 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */ 00027 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */ 00028 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */ 00029 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */ 00030 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */ 00031 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */ 00032 /* OTHER DEALINGS IN THE SOFTWARE. */ 00033 /* */ 00034 /************************************************************************/ 00035 00036 #ifndef VIGRA_RFFTW_HXX 00037 #define VIGRA_RFFTW_HXX 00038 00039 #include "array_vector.hxx" 00040 #include "fftw.hxx" 00041 #include <rfftw.h> 00042 00043 namespace vigra { 00044 00045 namespace detail { 00046 00047 struct FFTWSinCosConfig 00048 { 00049 ArrayVector<fftw_real> twiddles; 00050 ArrayVector<fftw_real> fftwInput; 00051 ArrayVector<fftw_complex> fftwTmpResult; 00052 fftw_real norm; 00053 rfftwnd_plan fftwPlan; 00054 }; 00055 00056 template <class SrcIterator, class SrcAccessor, 00057 class DestIterator, class DestAccessor, 00058 class Config> 00059 void 00060 cosineTransformLineImpl(SrcIterator s, SrcIterator send, SrcAccessor src, 00061 DestIterator d, DestAccessor dest, 00062 Config & config) 00063 { 00064 int size = send - s; 00065 int ns2 = size / 2; 00066 int nm1 = size - 1; 00067 int modn = size % 2; 00068 00069 if(size <= 0) 00070 return; 00071 00072 fftw_real const * twiddles = config.twiddles.begin(); 00073 fftw_real * fftwInput = config.fftwInput.begin(); 00074 fftw_complex * fftwTmpResult = config.fftwTmpResult.begin(); 00075 fftw_real norm = config.norm; 00076 rfftwnd_plan fftwPlan = config.fftwPlan; 00077 00078 switch(size) 00079 { 00080 case 1: 00081 { 00082 dest.set(src(s) / norm, d); 00083 break; 00084 } 00085 case 2: 00086 { 00087 dest.set((src(s) + src(s, 1)) / norm, d); 00088 dest.set((src(s) - src(s, 1)) / norm, d, 1); 00089 break; 00090 } 00091 case 3: 00092 { 00093 fftw_real x1p3 = src(s) + src(s, 2); 00094 fftw_real tx2 = 2.0 * src(s, 1); 00095 00096 dest.set((x1p3 + tx2) / norm, d, 0); 00097 dest.set((src(s) - src(s, 2)) / norm, d, 1); 00098 dest.set((x1p3 - tx2) / norm, d, 2); 00099 break; 00100 } 00101 default: 00102 { 00103 fftw_real c1 = src(s) - src(s, nm1); 00104 fftwInput[0] = src(s) + src(s, nm1); 00105 for(int k=1; k<ns2; ++k) 00106 { 00107 int kc = nm1 - k; 00108 fftw_real t1 = src(s, k) + src(s, kc); 00109 fftw_real t2 = src(s, k) - src(s, kc); 00110 c1 = c1 + twiddles[kc] * t2; 00111 t2 = twiddles[k] * t2; 00112 fftwInput[k] = t1 - t2; 00113 fftwInput[kc] = t1 + t2; 00114 } 00115 00116 if (modn != 0) 00117 { 00118 fftwInput[ns2] = 2.0*src(s, ns2); 00119 } 00120 rfftwnd_one_real_to_complex(fftwPlan, fftwInput, fftwTmpResult); 00121 dest.set(fftwTmpResult[0].re / norm, d, 0); 00122 for(int k=1; k<(size+1)/2; ++k) 00123 { 00124 dest.set(fftwTmpResult[k].re, d, 2*k-1); 00125 dest.set(fftwTmpResult[k].im, d, 2*k); 00126 } 00127 fftw_real xim2 = dest(d, 1); 00128 dest.set(c1 / norm, d, 1); 00129 for(int k=3; k<size; k+=2) 00130 { 00131 fftw_real xi = dest(d, k); 00132 dest.set(dest(d, k-2) - dest(d, k-1) / norm, d, k); 00133 dest.set(xim2 / norm, d, k-1); 00134 xim2 = xi; 00135 } 00136 if (modn != 0) 00137 { 00138 dest.set(xim2 / norm, d, size-1); 00139 } 00140 } 00141 } 00142 } 00143 00144 } // namespace detail 00145 00146 template <class SrcTraverser, class SrcAccessor, 00147 class DestTraverser, class DestAccessor> 00148 void cosineTransformX(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 00149 DestTraverser dul, DestAccessor dest, fftw_real norm) 00150 { 00151 int w = slr.x - sul.x; 00152 int h = slr.y - sul.y; 00153 00154 detail::FFTWSinCosConfig config; 00155 00156 // horizontal transformation 00157 int ns2 = w / 2; 00158 int nm1 = w - 1; 00159 int modn = w % 2; 00160 00161 config.twiddles.resize(w+1); 00162 config.fftwInput.resize(w+1); 00163 config.fftwTmpResult.resize(w+1); 00164 config.norm = norm; 00165 config.fftwPlan = rfftw2d_create_plan(1, nm1, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE ); 00166 00167 fftw_real dt = M_PI / nm1; 00168 for(int k=1; k<ns2; ++k) 00169 { 00170 fftw_real f = dt * k; 00171 config.twiddles[k] = 2.0*VIGRA_CSTD::sin(f); 00172 config.twiddles[nm1-k] = 2.0*VIGRA_CSTD::cos(f); 00173 } 00174 00175 for(; sul.y != slr.y; ++sul.y, ++dul.y) 00176 { 00177 typename SrcTraverser::row_iterator s = sul.rowIterator(); 00178 typename SrcTraverser::row_iterator send = s + w; 00179 typename DestTraverser::row_iterator d = dul.rowIterator(); 00180 cosineTransformLineImpl(s, send, src, d, dest, config); 00181 } 00182 00183 rfftwnd_destroy_plan(config.fftwPlan); 00184 } 00185 00186 template <class SrcTraverser, class SrcAccessor, 00187 class DestTraverser, class DestAccessor> 00188 void cosineTransformY(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 00189 DestTraverser dul, DestAccessor dest, fftw_real norm) 00190 { 00191 int w = slr.x - sul.x; 00192 int h = slr.y - sul.y; 00193 00194 detail::FFTWSinCosConfig config; 00195 00196 // horizontal transformation 00197 int ns2 = h / 2; 00198 int nm1 = h - 1; 00199 int modn = h % 2; 00200 00201 config.twiddles.resize(h + 1); 00202 config.fftwInput.resize(h + 1); 00203 config.fftwTmpResult.resize(h + 1); 00204 config.norm = norm; 00205 config.fftwPlan = rfftw2d_create_plan(1, nm1, FFTW_REAL_TO_COMPLEX, FFTW_ESTIMATE ); 00206 00207 fftw_real dt = M_PI / nm1; 00208 for(int k=1; k<ns2; ++k) 00209 { 00210 fftw_real f = dt * k; 00211 config.twiddles[k] = 2.0*VIGRA_CSTD::sin(f); 00212 config.twiddles[nm1-k] = 2.0*VIGRA_CSTD::cos(f); 00213 } 00214 00215 for(; sul.x != slr.x; ++sul.x, ++dul.x) 00216 { 00217 typename SrcTraverser::column_iterator s = sul.columnIterator(); 00218 typename SrcTraverser::column_iterator send = s + h; 00219 typename DestTraverser::column_iterator d = dul.columnIterator(); 00220 cosineTransformLineImpl(s, send, src, d, dest, config); 00221 } 00222 00223 rfftwnd_destroy_plan(config.fftwPlan); 00224 } 00225 00226 template <class SrcTraverser, class SrcAccessor, 00227 class DestTraverser, class DestAccessor> 00228 inline void 00229 realFourierTransformXEvenYEven(SrcTraverser sul, SrcTraverser slr, SrcAccessor src, 00230 DestTraverser dul, DestAccessor dest, fftw_real norm) 00231 { 00232 BasicImage<fftw_real> tmp(slr - sul); 00233 cosineTransformX(sul, slr, src, tmp.upperLeft(), tmp.accessor(), 1.0); 00234 cosineTransformY(tmp.upperLeft(), tmp.lowerRight(), tmp.accessor(), dul, dest, norm); 00235 } 00236 00237 template <class SrcTraverser, class SrcAccessor, 00238 class DestTraverser, class DestAccessor> 00239 inline void 00240 realFourierTransformXEvenYEven(triple<SrcTraverser, SrcTraverser, SrcAccessor> src, 00241 pair<DestTraverser, DestAccessor> dest, fftw_real norm) 00242 { 00243 realFourierTransformXEvenYEven(src.first, src.second, src.third, dest.first, dest.second, norm); 00244 } 00245 00246 } // namespace vigra 00247 00248 #endif // VIGRA_RFFTW_HXX
© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de) |
html generated using doxygen and Python
|