1

Basically I am trying to implement a very basic version of the Wiener filter on a grey scale image, using the a stripped down Wiener equation: (1/(SNR))*DFT(Image) after which I take the IDFT of the whole thing. My problem is that my output image which is supposed to be filtered looks exactly like the input, and therefore it seems that the pixel values aren't changing at all. Can anyone please indicate to me where I am going wrong? Here's the code I am currently using:

#include "opencv2/imgproc/imgproc.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv/cv.hpp"
#include "opencv/cxcore.hpp"
#include <iostream>


using namespace cv;
using namespace std;

void updateMag(Mat complex);
Mat updateResult(Mat complex);

Mat computeDFT(Mat image);
Mat DFT2(Mat I);
void shift(Mat magI);

int kernel_size = 0;

int main( int argc, char** argv )
{

    Mat  result;
    String file;
    file = " << SAMPLE FILE >>";

    Mat image = imread("/Users/John/Desktop/house.png", CV_LOAD_IMAGE_GRAYSCALE);
    namedWindow( "Orginal window", CV_WINDOW_AUTOSIZE  );// Create a window for display.
    imshow( "Orginal window", image );                   // Show our image inside it.

    float x = 1/0.001;
    Mat complex = computeDFT(image); // DFT of image

    updateMag(complex);         // compute magnitude of complex, switch to logarithmic scale and display...

    Mat fourierImage(complex.size(), complex.type());
    fourierImage = cv::Scalar::all(x);


    //cout<< "Fourier = " << endl << fourierImage << endl;
    //Mat complexFourier = computeDFT(fourierImage);
    //cout << "1" << endl << complexFourier.type() << endl << complexFourier.type() << endl;

    //complex = complex.mul(fourierImage);
    //mulSpectrums(complex, fourierImage, complex, DFT_ROWS);

    complex = complex.mul(x);
    result = updateResult(complex);      // do inverse transform and display the result image

    waitKey(0);

    return 0;
}



Mat updateResult(Mat complex)
{
    Mat work;
    //work.convertTo(work, CV_32F);
    idft(complex, work);
    //dft(complex, work, DFT_INVERSE + DFT_SCALE);
    Mat planes[] = {Mat::zeros(complex.size(), complex.type()), Mat::zeros(complex.size(), complex.type())};
    split(work, planes);                // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))

    magnitude(planes[0], planes[1], work);    // === sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)
    normalize(work, work, 1, 0, NORM_MINMAX);
    imshow("result", work);
    return work;
}

void updateMag(Mat complex )
{

    Mat magI;
    Mat planes[] = {Mat::zeros(complex.size(), CV_32F), Mat::zeros(complex.size(), CV_32F)};
    split(complex, planes);                // planes[0] = Re(DFT(I)), planes[1] = Im(DFT(I))

    magnitude(planes[0], planes[1], magI);    // sqrt(Re(DFT(I))^2 + Im(DFT(I))^2)

    // switch to logarithmic scale: log(1 + magnitude)
    magI += Scalar::all(1);
    log(magI, magI);

    shift(magI);
    normalize(magI, magI, 1, 0, NORM_INF); // Transform the matrix with float values into a
    // viewable image form (float between values 0 and 1).
    imshow("spectrum", magI);
}


Mat computeDFT(Mat image) {
    Mat padded;                            //expand input image to optimal size
    int m = getOptimalDFTSize( image.rows );
    int n = getOptimalDFTSize( image.cols ); // on the border add zero values
    copyMakeBorder(image, padded, 0, m - image.rows, 0, n - image.cols, BORDER_CONSTANT, Scalar::all(0));
    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complex;
    merge(planes, 2, complex);         // Add to the expanded another plane with zeros
    dft(complex, complex, DFT_COMPLEX_OUTPUT);  // furier transform
    return complex;
}



void shift(Mat magI) {

    // crop if it has an odd number of rows or columns
    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

    int cx = magI.cols/2;
    int cy = magI.rows/2;

    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

    Mat tmp;                            // swap quadrants (Top-Left with Bottom-Right)
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);                     // swap quadrant (Top-Right with Bottom-Left)
    q2.copyTo(q1);
    tmp.copyTo(q2);
}

Mat DFT2(Mat I)
{
    Mat padded;                            //expand input image to optimal size
    int m = getOptimalDFTSize( I.rows );
    int n = getOptimalDFTSize( I.cols ); // on the border add zero values
    copyMakeBorder(I, padded, 0, m - I.rows, 0, n - I.cols, BORDER_CONSTANT, Scalar::all(0));

    Mat planes[] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)};
    Mat complexI;
    merge(planes, 2, complexI);         // Add to the expanded another plane with zeros

    dft(complexI, complexI);            // this way the result may fit in the source matrix

    // compute the magnitude and switch to logarithmic scale
    // => log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    split(complexI, planes);                   // planes[0] = Re(DFT(I), planes[1] = Im(DFT(I))
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude
    Mat magI = planes[0];

    magI += Scalar::all(1);                    // switch to logarithmic scale
    log(magI, magI);

    // crop the spectrum, if it has an odd number of rows or columns
    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));

    // rearrange the quadrants of Fourier image  so that the origin is at the image center
    int cx = magI.cols/2;
    int cy = magI.rows/2;

    Mat q0(magI, Rect(0, 0, cx, cy));   // Top-Left - Create a ROI per quadrant
    Mat q1(magI, Rect(cx, 0, cx, cy));  // Top-Right
    Mat q2(magI, Rect(0, cy, cx, cy));  // Bottom-Left
    Mat q3(magI, Rect(cx, cy, cx, cy)); // Bottom-Right

    Mat tmp;                           // swap quadrants (Top-Left with Bottom-Right)
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);

    q1.copyTo(tmp);                    // swap quadrant (Top-Right with Bottom-Left)
    q2.copyTo(q1);
    tmp.copyTo(q2);

    normalize(magI, magI, 0, 1, CV_MINMAX); // Transform the matrix with float values into a
    // viewable image form (float between values 0 and 1).


    return complexI;
}
  • If you are interested in working with python, [HAVE A LOOK HERE](http://scikit-image.org/docs/dev/api/skimage.filters.html#skimage.filters.wiener) – Jeru Luke Feb 23 '17 at 17:34
  • @JeruLuke Thanks for the comment however I don't really have a choice as this is for a thesis and I have to use C++ – John Napier Feb 23 '17 at 22:30
  • hmmm... don't have hopes on me... but I'll try my best. I am interested in this solution too – Jeru Luke Feb 24 '17 at 03:58
  • 1
    @JeruLuke I have been looking around and i found this post very useful - Here is the [link](http://stackoverflow.com/questions/10269456/inverse-fourier-transformation-in-opencv) – John Napier Feb 27 '17 at 10:55
  • Thnx for letting me know :) – Jeru Luke Feb 27 '17 at 11:55

0 Answers0