Converting fftw library calls to OpenCV calls

Question

I am trying to use OpenCV instead of fftw because of the more permissive license of OpenCV. With the help of FFTW fftwf_plan_r2r_2d() with FFTW_REDFT01 equivalent

I am aware that it is currently wrong, it's just what I have as an initial idea. I need some guidance to get exact functionality of fftw with OpenCV, I am offering bounty for someone who can help me do that.

I only need those particular functions and I am not interested in converting other fftw functions for now.

Thanks!

Anyone?

h:

#pragma once

#define DO_NOT_USE_FFTW

#ifdef DO_NOT_USE_FFTW
enum fftwf_r2r_kind
{
    FFTW_R2HC,    //computes a real-input DFT with output in “halfcomplex” format, i.e. real and imaginary parts for a transform of size n stored as: r0, r1, r2, ..., rn/2, i(n+1)/2-1, ..., i2, i1 (Logical N=n, inverse is FFTW_HC2R.)
    FFTW_HC2R,    //computes the reverse of FFTW_R2HC, above. (Logical N=n, inverse is FFTW_R2HC.)
    FFTW_DHT,     //computes a discrete Hartley transform. (Logical N=n, inverse is FFTW_DHT.)
    FFTW_REDFT00, //computes an REDFT00 transform, i.e. a DCT-I. (Logical N=2*(n-1), inverse is FFTW_REDFT00.)
    FFTW_REDFT10, //computes an REDFT10 transform, i.e. a DCT-II (sometimes called “the” DCT). (Logical N=2*n, inverse is FFTW_REDFT01.)
    FFTW_REDFT01, //computes an REDFT01 transform, i.e. a DCT-III (sometimes called “the” IDCT, being the inverse of DCT-II). (Logical N=2*n, inverse is FFTW_REDFT=10.)
    FFTW_REDFT11, //computes an REDFT11 transform, i.e. a DCT-IV. (Logical N=2*n, inverse is FFTW_REDFT11.)
    FFTW_RODFT00, //computes an RODFT00 transform, i.e. a DST-I. (Logical N=2*(n+1), inverse is FFTW_RODFT00.)
    FFTW_RODFT10, //computes an RODFT10 transform, i.e. a DST-II. (Logical N=2*n, inverse is FFTW_RODFT01.)
    FFTW_RODFT01, //computes an RODFT01 transform, i.e. a DST-III. (Logical N=2*n, inverse is FFTW_RODFT=10.)
    FFTW_RODFT11, //computes an RODFT11 transform, i.e. a DST-IV. (Logical N=2*n, inverse is FFTW_RODFT11.)
};

struct fftwf_plan_imp
{
     int nx;
     int ny;
     float * in;
     float * out; 
     fftwf_r2r_kind kindx;
     fftwf_r2r_kind kindy;
     unsigned flags;
};

#define fftwf_plan         fftwf_plan_imp *
#define fftwf_malloc(x)    malloc(x)
#define fftwf_free(x)      free(x)
#define FFTW_DESTROY_INPUT 1
#define FFTW_ESTIMATE      1

fftwf_plan fftwf_plan_r2r_2d(int nx, int ny, float * in, float * out, fftwf_r2r_kind kindx, fftwf_r2r_kind kindy, unsigned flags);
void fftwf_execute(fftwf_plan x);
void fftwf_destroy_plan(fftwf_plan x);
void fftwf_cleanup();
#else
#include <fftw3.h>
#endif

src:

#include "KissInterface.h"
#include "kiss_fftr.h"
#include "kiss_fftndr.h"
#include "opencv.hpp"

//currently cv::dct supports even-size arrays (2, 4, 6 ...). For data analysis and approximation you can pad the array when necessary.

#ifdef DO_NOT_USE_FFTW
fftwf_plan fftwf_plan_r2r_2d(int nx, int ny, float * in, float * out, fftwf_r2r_kind kindx, fftwf_r2r_kind kindy, unsigned flags)
{
    fftwf_plan plan = new fftwf_plan_imp;

    plan->nx    = nx;
    plan->ny    = ny;
    plan->in    = in;
    plan->out   = out;
    plan->kindx = kindx;
    plan->kindy = kindy;
    plan->flags = flags;

    return plan;
}

void fftwf_execute(fftwf_plan x)
{
    int newx = x->nx;
    int newy = x->ny;
    if(newx % 2 != 0)++newx;
    if(newy % 2 != 0)++newy;

    cv::Mat Input (cv::Size(newx, newy), CV_32FC1, x->in); //AUTO_STEP or not?
    cv::Mat Output(cv::Size(newx, newy), CV_32FC1, x->out); //AUTO_STEP or not?

    if(x->kindx == FFTW_REDFT00 && x->kindy == FFTW_REDFT00)
    {
        cv::dct(Input, Output);
    }
    else if(x->kindx == FFTW_REDFT10 && x->kindy == FFTW_REDFT10)
    {
        cv::dct(Input, Output);

        Output *= (4 * sqrt(float(newx / 2)) * sqrt(float(newy / 2)));
        Output.row(0) *= 1.41421356237;
        Output.col(0) *= 1.41421356237;
    }
    else if(x->kindx == FFTW_REDFT01 && x->kindy == FFTW_REDFT01)
    {
        // First re-scale the data for idct():
        Input /= (4 * sqrt(float(newx / 2)) * sqrt(float(newy / 2)));
        Input.row(0) /= 1.41421356237;
        Input.col(0) /= 1.41421356237;

        cv::idct(Input, Output); // this will return the input exactly

        // However, the transforms computed by FFTW are unnormalized, exactly like the corresponding, 
        // so computing a transform followed by its inverse yields the original array scaled by N, where N is the logical DFT size. 
        // The logical DFT size: Logical N=2*n for each axis, this is th implicit symmetrization
        // of the image: reflect right and then reflect both halves down.
        int logicalSizeN = (2 * newx) * (2 * newy);
        Output *= logicalSizeN; // scale to be like FFTW result
    }

    memcpy(x->out, Output.ptr<float>(0), x->nx * x->ny);
}

void fftwf_destroy_plan(fftwf_plan x)
{
    delete x;
    x = NULL;
}

void fftwf_cleanup()
{

}
#endif 

Answer 1

OpenCV's FFT/DCT works like Matlab's implementations. Knowing this may help you find better references and answers.
In my question and answer you have several links that may help get you started.

Recently I had to go back to this code and rework it and it turned out that since I wanted to use the transforms and then transform back, I could just avoid all the extra scaling, before and after, and it just worked.
Of course, this means that the intermediate values (between the transforms) are not the same as with fftw, but the results after the inverse-transform were.