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