[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | 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 /*    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)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
VIGRA 1.6.0 (5 Nov 2009)