Обратное преобразование Фурье ручное слияние реального и воображаемого

Question

Я создаю пару методов-оберток для преобразования изображения с четырьмя изображениями и столкнулся с проблемой. Всякий раз, когда я использую стандартные методы:

void fft::dfft(Mat* img,Mat* result){ 
    Mat image; 
    img->convertTo(image, CV_32F); 
    cv::dft(image, *result, cv::DFT_SCALE|cv::DFT_COMPLEX_OUTPUT); 
} 

Обратить операции:

void fft::idft(Mat* fourier, Mat* img) { 
    //invert: 
    Mat inverseTransform; 
    dft(*fourier, inverseTransform, cv::DFT_INVERSE|cv::DFT_REAL_OUTPUT); 
    //restore image to 8 bits (fourier needs to use float at least 
    Mat res; 
    inverseTransform.convertTo(res, CV_8U); 
    res.copyTo(*img); 
} 

все в порядке. Однако, когда я пытаюсь использовать эту пару:

void fft::dfft(Mat* img, Mat* mag, Mat* phi, Mat* re, Mat* im) { 

/*SPECAL NOTE: dft is faster for images with size being multiple of two three and five. An OCV fmethod getOptrimalDFTSize comes in handy*/ 
Mat paddedImg; //expand input image to the optimal size 
int m = getOptimalDFTSize(img->rows); 
int n = getOptimalDFTSize(img->cols); 
copyMakeBorder(*img, paddedImg, 0, m - img->rows, 0, n-img->cols, BORDER_CONSTANT, Scalar::all(0)); 

/*prepare space to store the image*/ 
Mat planes[] = {Mat_<float>(paddedImg), Mat::zeros(paddedImg.size(), CV_32F)}; 
Mat complexImg; 
merge(planes, 2, complexImg); 

/*Actual dft:*/ 
dft(complexImg, complexImg); 

/*Get magnitude:*/ 
split(complexImg, planes); 
Mat magImg, phasImg; 
planes[0].copyTo(*re); 
planes[1].copyTo(*im); 

cout << "intern Real = "<< planes[0] << endl; 

magnitude(planes[0], planes[1], magImg); //magnitude will be stored in planes[0] 
phase(planes[0], planes[1], phasImg); 


magImg.copyTo(*mag); 
phasImg.copyTo(*phi); 
#ifdef debug 
namedWindow("Input Image", 0); 
imshow("Input Image", *img); 
namedWindow("Image magnitude", 0); 
imshow("Image magnitude", magImg); 
waitKey(0); 
#endif 
} 

Чтобы извлечь конкретные параметры, я не могу вернуть их обратно. В настоящее время я сосредотачиваюсь на Real и Imaginary part, потому что я чувствую, что Mag/Phi будет работать одинаково, только сначала преобразуясь в реальную и мнимую назад из Mag/Phi.

void fft::idft(Mat* re, Mat* im, Mat* img, bool dummy) { 
Mat inverseTransform; 

Mat fourier(re->rows, re->cols, CV_32FC2); 
vector<Mat> planes; 
planes.push_back(*re); 
planes.push_back(*im); 

// повторно и им попасть внутрь правильно (проверено)

Mat padded; 
Mat complex; 
Mat pp[] = {Mat_<float>(*re), Mat::zeros(re->size(), CV_32F)}; 
pp[0] = *re; 
pp[1] = *im; 
merge(pp, 2, complex);   // Add to the expanded another plane with zeros 

cv::merge(planes, fourier); 

dft(complex, inverseTransform, cv::DFT_INVERSE|cv::DFT_REAL_OUTPUT); 

Mat res; 
inverseTransform.convertTo(res, CV_8U); 
cv::imshow("iFFT1", res); 
dft(fourier, inverseTransform, cv::DFT_INVERSE|cv::DFT_REAL_OUTPUT); 
inverseTransform.convertTo(res, CV_8U); 
cv::imshow("iFFT2", res); 


res.copyTo(*img); 
} 

Оба подхода (с вектором и матрицы С.) производят один и тот же белый прямоугольник. Чтобы узнать, не делает ли это что-либо, я прошел fft :: idft (& re, & re, & img, true); В этом случае результат был идентичным в обоих окнах (некоторый случайный мусор, но идентичный)

Я подозреваю, что проблема заключается в неправильной процедуре слияния, но у меня нет других идей по восстановлению исходной формы?

Заранее благодарим за понимание.

Answer 1

Вот мой фильтр домен Фурье, я думаю, что это должно быть полезно:

#pragma once 
#include <string> 
#include <iostream> 
#include "opencv2/opencv.hpp" 
using namespace std; 
using namespace cv; 
//---------------------------------------------------------- 
// Recombinate quadrants 
//---------------------------------------------------------- 
void Recomb(Mat &src,Mat &dst) 
{ 
    int cx = src.cols>>1; 
    int cy = src.rows>>1; 
    Mat tmp; 
    tmp.create(src.size(),src.type()); 
    src(Rect(0, 0, cx, cy)).copyTo(tmp(Rect(cx, cy, cx, cy))); 
    src(Rect(cx, cy, cx, cy)).copyTo(tmp(Rect(0, 0, cx, cy))); 
    src(Rect(cx, 0, cx, cy)).copyTo(tmp(Rect(0, cy, cx, cy))); 
    src(Rect(0, cy, cx, cy)).copyTo(tmp(Rect(cx, 0, cx, cy))); 
    dst=tmp; 
} 
//---------------------------------------------------------- 
// Forward fft 
//---------------------------------------------------------- 
void ForwardFFT(Mat &Src, Mat *FImg) 
{ 
    int M = getOptimalDFTSize(Src.rows); 
    int N = getOptimalDFTSize(Src.cols); 
    Mat padded;  
    copyMakeBorder(Src, padded, 0, M - Src.rows, 0, N - Src.cols, BORDER_CONSTANT, Scalar::all(0)); 
    // Create complex image 
    // planes[0] image , planes[1] filled by zeroes 
    Mat planes[2] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}; 
    Mat complexImg; 
    merge(planes, 2, complexImg); 
    dft(complexImg, complexImg,DFT_SCALE);  
    // After tranform we also get complex image 
    split(complexImg, planes); 

    // 
    planes[0] = planes[0](Rect(0, 0, planes[0].cols & -2, planes[0].rows & -2)); 
    planes[1] = planes[1](Rect(0, 0, planes[1].cols & -2, planes[1].rows & -2)); 

    Recomb(planes[0],planes[0]); 
    Recomb(planes[1],planes[1]); 

    FImg[0]=planes[0].clone(); 
    FImg[1]=planes[1].clone(); 
} 
//---------------------------------------------------------- 
// Inverse FFT 
//---------------------------------------------------------- 
void InverseFFT(Mat *FImg,Mat &Dst) 
{ 
    Recomb(FImg[0],FImg[0]); 
    Recomb(FImg[1],FImg[1]); 
    Mat complexImg; 
    merge(FImg, 2, complexImg); 
    // Inverse transform 
    dft(complexImg, complexImg, DFT_INVERSE); 
    split(complexImg, FImg);   
    FImg[0].copyTo(Dst); 
} 
//---------------------------------------------------------- 
// Forward FFT using Magnitude and phase 
//---------------------------------------------------------- 
void ForwardFFT_Mag_Phase(Mat &src, Mat &Mag,Mat &Phase) 
{ 
    Mat planes[2]; 
    ForwardFFT(src,planes); 
    Mag.zeros(planes[0].rows,planes[0].cols,CV_32F); 
    Phase.zeros(planes[0].rows,planes[0].cols,CV_32F); 
    cv::cartToPolar(planes[0],planes[1],Mag,Phase); 
} 
//---------------------------------------------------------- 
// Inverse FFT using Magnitude and phase 
//---------------------------------------------------------- 
void InverseFFT_Mag_Phase(Mat &Mag,Mat &Phase, Mat &dst) 
{ 
    Mat planes[2]; 
    planes[0].create(Mag.rows,Mag.cols,CV_32F); 
    planes[1].create(Mag.rows,Mag.cols,CV_32F); 
    cv::polarToCart(Mag,Phase,planes[0],planes[1]); 
    InverseFFT(planes,dst); 
} 
//---------------------------------------------------------- 
// MAIN 
//---------------------------------------------------------- 
int main(int argc, char* argv[]) 
{ 
    // src image 
    Mat img; 
    // Magnitude 
    Mat Mag; 
    // Phase 
    Mat Phase; 
    // Image loading 
    img=imread("d:\\ImagesForTest\\lena.jpg",0); 
    resize(img,img,Size(512,512)); 
    // Image size 
    cout<<img.size().width<<endl; 
    cout<<img.size().height<<endl; 

    // 
    ForwardFFT_Mag_Phase(img,Mag,Phase);  

    //---------------------------------------------------------- 
    // Filter 
    //---------------------------------------------------------- 
    // draw ring  
    int R=100; // External radius 
    int r=30; // internal radius 
    Mat mask; 
    mask.create(Mag.cols,Mag.rows,CV_32F); 
    int cx = Mag.cols>>1; 
    int cy = Mag.rows>>1;  
    mask=1; 
    cv::circle(mask,cv::Point(cx,cy),R,CV_RGB(0,0,0),-1); 
    cv::circle(mask,cv::Point(cx,cy),r,CV_RGB(1,1,1),-1); 
    //mask=1-mask; // uncomment for inversion 

    //cv::multiply(Mag,mask,Mag); // uncomment to turn filter on 
    //cv::multiply(Phase,mask,Phase); 
    //---------------------------------------------------------- 
    // Inverse transform 
    //---------------------------------------------------------- 
    InverseFFT_Mag_Phase(Mag,Phase,img);  

    //---------------------------------------------------------- 
    // Results output 
    //---------------------------------------------------------- 
    // 
    Mat LogMag; 
    LogMag.zeros(Mag.rows,Mag.cols,CV_32F); 
    LogMag=(Mag+1); 
    cv::log(LogMag,LogMag); 
    //--------------------------------------------------- 
    imshow("Magnitude Log", LogMag); 
    imshow("Phase", Phase); 
    // img - now in CV_32FC1 format,we need CV_8UC1 or scale it by factor 1.0/255.0 
    img.convertTo(img,CV_8UC1); 
    imshow("Filtering result", img);  
    //---------------------------------------------------------- 
    // Wait key press 
    //---------------------------------------------------------- 
    waitKey(0); 
    return 0; 
    } 

К сожалению, я не могу воспроизвести вам результат, я проверил код и он работает отлично:

#pragma once 
#include <string> 
#include <iostream> 
#include "opencv2/opencv.hpp" 
using namespace std; 
using namespace cv; 
//---------------------------------------------------------- 
// Recombinate quadrants 
//---------------------------------------------------------- 
void Recomb(Mat &src,Mat &dst) 
{ 
    int cx = src.cols>>1; 
    int cy = src.rows>>1; 
    Mat tmp; 
    tmp.create(src.size(),src.type()); 
    src(Rect(0, 0, cx, cy)).copyTo(tmp(Rect(cx, cy, cx, cy))); 
    src(Rect(cx, cy, cx, cy)).copyTo(tmp(Rect(0, 0, cx, cy))); 
    src(Rect(cx, 0, cx, cy)).copyTo(tmp(Rect(0, cy, cx, cy))); 
    src(Rect(0, cy, cx, cy)).copyTo(tmp(Rect(cx, 0, cx, cy))); 
    dst=tmp; 
} 
//---------------------------------------------------------- 
// Forward fft 
//---------------------------------------------------------- 
void ForwardFFT(Mat &Src, Mat *FImg) 
{ 
    int M = getOptimalDFTSize(Src.rows); 
    int N = getOptimalDFTSize(Src.cols); 
    Mat padded;  
    copyMakeBorder(Src, padded, 0, M - Src.rows, 0, N - Src.cols, BORDER_CONSTANT, Scalar::all(0)); 
    // Create complex image 
    // planes[0] image , planes[1] filled by zeroes 
    Mat planes[2] = {Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F)}; 
    Mat complexImg; 
    merge(planes, 2, complexImg); 
    dft(complexImg, complexImg,DFT_SCALE);  
    // After tranform we also get complex image 
    split(complexImg, planes); 

    // 
    planes[0] = planes[0](Rect(0, 0, planes[0].cols & -2, planes[0].rows & -2)); 
    planes[1] = planes[1](Rect(0, 0, planes[1].cols & -2, planes[1].rows & -2)); 

    Recomb(planes[0],planes[0]); 
    Recomb(planes[1],planes[1]); 

    FImg[0]=planes[0].clone(); 
    FImg[1]=planes[1].clone(); 
} 
//---------------------------------------------------------- 
// Inverse FFT 
//---------------------------------------------------------- 
void InverseFFT(Mat *FImg,Mat &Dst) 
{ 
    Recomb(FImg[0],FImg[0]); 
    Recomb(FImg[1],FImg[1]); 
    Mat complexImg; 
    merge(FImg, 2, complexImg); 
    // Inverse transform 
    dft(complexImg, complexImg, DFT_INVERSE); 
    split(complexImg, FImg);   
    FImg[0].copyTo(Dst); 
} 
//---------------------------------------------------------- 
// Forward FFT using Magnitude and phase 
//---------------------------------------------------------- 
void ForwardFFT_Mag_Phase(Mat &src, Mat &Mag,Mat &Phase) 
{ 
    Mat planes[2]; 
    ForwardFFT(src,planes); 
    Mag.zeros(planes[0].rows,planes[0].cols,CV_32F); 
    Phase.zeros(planes[0].rows,planes[0].cols,CV_32F); 
    cv::cartToPolar(planes[0],planes[1],Mag,Phase); 
} 
//---------------------------------------------------------- 
// Inverse FFT using Magnitude and phase 
//---------------------------------------------------------- 
void InverseFFT_Mag_Phase(Mat &Mag,Mat &Phase, Mat &dst) 
{ 
    Mat planes[2]; 
    planes[0].create(Mag.rows,Mag.cols,CV_32F); 
    planes[1].create(Mag.rows,Mag.cols,CV_32F); 
    cv::polarToCart(Mag,Phase,planes[0],planes[1]); 
    InverseFFT(planes,dst); 
} 
//---------------------------------------------------------- 
// MAIN 
//---------------------------------------------------------- 
int main(int argc, char* argv[]) 
{ 
    // src image 
    Mat img; 
    // Magnitude 
    Mat Mag; 
    // Phase 
    Mat Phase; 
    // Image loading (grayscale) 
    img=imread("d:\\ImagesForTest\\lena.jpg",0); 
    ForwardFFT_Mag_Phase(img,Mag,Phase);  
    //---------------------------------------------------------- 
    // Inverse transform 
    //---------------------------------------------------------- 
    InverseFFT_Mag_Phase(Mag,Phase,img);  
    img.convertTo(img,CV_8UC1); 
    imshow("Filtering result", img);  
    //---------------------------------------------------------- 
    // Wait key press 
    //---------------------------------------------------------- 
    waitKey(0); 
    return 0; 
} 

И он получает то же изображение, что и вход.