【图像处理笔记】傅里叶变换

#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

int main() {    
    clock_t start = clock();
    Mat src = imread("E:/C++projects/CH11/CH11/5.png", 0);
    imshow("src", src);
    // 填充到傅里叶变换最佳尺寸 （gpu下加速明显）
    int h = getOptimalDFTSize(src.rows);
    int w = getOptimalDFTSize(src.cols);
    Mat padded;
    copyMakeBorder(src, padded, 0, h - src.rows, 0, w - src.cols, BORDER_CONSTANT, Scalar::all(0));
    Mat pMat = Mat(padded.size(), CV_32FC1);
    
    // 中心化 （在时域做中心化，傅里叶变换后的频谱图低频在中心）
    for (int i = 0; i < padded.rows; i++)
    {
        for (int j = 0; j < padded.cols; j++)
        {
            pMat.at<float>(i, j) = pow(-1, i + j) * padded.at<uchar>(i, j);
        }
    }    
    
    // 傅里叶变换
    Mat planes[] = { pMat, Mat::zeros(padded.size(), CV_32F) };
    Mat complexImg;
    merge(planes, 2, complexImg);
    dft(complexImg, complexImg, DFT_COMPLEX_OUTPUT);

    // 显示频谱图
    split(complexImg, planes);
    Mat magMat;
    magnitude(planes[0], planes[1], magMat);
    magMat += 1;
    log(magMat, magMat);
    normalize(magMat, magMat, 0, 255, NORM_MINMAX);
    magMat.convertTo(magMat, CV_8UC1);

    // 低通滤波
    //Mat mask(Size(w, h), CV_32FC2, Scalar(0, 0));
    //circle(mask, Point(w / 2, h / 2), 50, Scalar(255, 255), -1);
    //complexImg = complexImg.mul(mask);

    //高通滤波
    Mat mask(Size(w, h), CV_32FC2, Scalar(255, 255));
    circle(mask, Point(w / 2, h / 2), 50, Scalar(0, 0), -1);
    complexImg = complexImg.mul(mask);

    // 傅里叶反变换
    Mat dst, tmp;
    idft(complexImg, tmp, DFT_REAL_OUTPUT);
    
    // 去中心化
    for (int i = 0; i < padded.rows; i++)
    {
        for (int j = 0; j < padded.cols; j++)
        {
            pMat.at<float>(i, j) = pow(-1, i + j) * tmp.at<float>(i, j);
        }
    }
    dst = pMat({ 0,0,src.cols ,src.rows }).clone();
    normalize(dst, dst, 0, 1, NORM_MINMAX);
    dst.convertTo(dst, CV_8U, 255);
    imshow("dst", dst);
    clock_t end = clock(); 
    cout << "spend time:" << end - start << "ms" << endl;
    waitKey(0);
    return 0; 
}

#include <opencv2/opencv.hpp>
using namespace cv;
using namespace std;

double point2Line(Point2f pointP, Point2f pointA, Point2f pointB)
{
    //求直线方程
    double A = 0, B = 0, C = 0;
    A = pointA.y - pointB.y;
    B = pointB.x - pointA.x;
    C = pointA.x * pointB.y - pointA.y * pointB.x;
    //代入点到直线距离公式
    double distance = 0;
    double tmp1 = abs(A * pointP.x + B * pointP.y + C);
    double tmp2 = ((float)sqrtf(A * A + B * B));
    distance = tmp1 / tmp2;
    return distance;
}
int main() {
    vector<string> filenames;
    glob("./", filenames);
    for (int n = 0; n < filenames.size(); n++) {
        string filename = filenames[n];
        Mat src = imread(filename, 0);
        // 填充到傅里叶变换的最佳尺寸
        int h = getOptimalDFTSize(src.rows);
        int w = getOptimalDFTSize(src.cols);
        Mat padded;
        copyMakeBorder(src, padded, 0, h - src.rows, 0, w - src.cols, BORDER_CONSTANT, Scalar::all(0));
        // 傅里叶变换
        Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
        Mat complexImg;
        merge(planes, 2, complexImg);
        dft(complexImg, complexImg, DFT_COMPLEX_INPUT | DFT_COMPLEX_OUTPUT);
        // 求幅值
        split(complexImg, planes);
        Mat magMat;
        magnitude(planes[0], planes[1], magMat);
        magMat += 1;
        log(magMat, magMat);
        normalize(magMat, magMat, 0, 255, NORM_MINMAX);
        magMat.convertTo(magMat, CV_8UC1);
        // 把零频移到中心
        Mat magImg = magMat.clone();
        int cx = magImg.cols / 2;
        int cy = magImg.rows / 2;
        Mat q1 = magImg({ 0, 0, cx, cy });
        Mat q2 = magImg({ 0, cy, cx, cy });
        Mat q3 = magImg({ cx, 0, cx, cy });
        Mat q4 = magImg({ cx, cy, cx, cy });
        Mat temp;
        q1.copyTo(temp);
        q4.copyTo(q1);
        temp.copyTo(q4);
        q2.copyTo(temp);
        q3.copyTo(q2);
        temp.copyTo(q3);
        // 霍夫直线检测有角度的斜线
        Mat binImg;
        threshold(magImg, binImg, 150, 255, THRESH_BINARY);
        Mat markImg;
        cvtColor(binImg, markImg, COLOR_GRAY2BGR);
        vector<Vec4i> lines;
        Vec4i l;
        HoughLinesP(binImg, lines, 1, CV_PI / 180.0, 30, 200, 50);
        Point2f p = Point2f(magImg.cols / 2.0, magImg.rows / 2.0);
        for (int i = 0; i < lines.size(); i++)
        {
            if (abs(lines[i][1] - lines[i][3]) > 15) {
                line(markImg, Point(lines[i][0], lines[i][1]), Point(lines[i][2], lines[i][3]), Scalar(0, 255, 0), 1, 8, 0);
                l = lines[i];
            }
        }
        float theta = atan((l[1] - l[3]) * 1.0 / (l[0] - l[2]) * src.cols / src.rows) * 180 / CV_PI;
        float angle = theta <= 0 ? theta + 90 : theta - 90;
        // 放射变换
        Point2f center = Point2f(src.cols / 2.0, src.rows / 2.0);
        Mat rotMat = getRotationMatrix2D(center, angle, 1.0);
        Mat dst = Mat::ones(src.size(), CV_8UC1);
        warpAffine(src, dst, rotMat, src.size(), 1, 0, Scalar(255, 255, 255));
        imshow("src", src);
        imshow("markImg", markImg);
        imshow("dst", dst);
        waitKey(0);
    }
    return 0;
}

...
Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
Mat complexImg;
merge(planes, 2, complexImg);
dft(complexImg, complexImg, DFT_COMPLEX_OUTPUT);
split(complexImg, planes);
Mat magMat;
magnitude(planes[0], planes[1], magMat);
// 相角重建
//divide(planes[0], magMat, planes[0]);
//divide(planes[1], magMat, planes[1]);
//谱重建
planes[0] = magMat;
planes[1] = Mat::zeros(padded.size(), CV_32F);
merge(planes, 2, complexImg);
Mat dst, tmp;
idft(complexImg, tmp, DFT_REAL_OUTPUT);
...