Преобразуване на извиквания на библиотеката fftw в извиквания на OpenCV

Опитвам се да използвам OpenCV вместо fftw поради по-разрешителния лиценз на OpenCV. С помощта на FFTW fftwf_plan_r2r_2d() с FFTW_REDFT01 еквивалент

Наясно съм, че в момента не е наред, просто това е, което имам като първоначална идея. Имам нужда от някои насоки, за да получа точната функционалност на fftw с OpenCV, предлагам награда за някой, който може да ми помогне да направя това.

Имам нужда само от тези конкретни функции и засега не се интересувам от конвертиране на други функции на fftw.

Благодаря!

Някой?

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