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