[ 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 /*    ( Version 1.6.0, Aug 13 2008 )                                    */
    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 /*        ullrich.koethe@iwr.uni-heidelberg.de    or                    */
    00012 /*        vigra@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 (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 (13 Aug 2008)