[ 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) Cognitive Systems Group, University of Hamburg, Germany | html generated using doxygen and Python VIGRA 1.3.2 (27 Jan 2005)