低通滤波

#pragma 
#include <opencv2/opencv.hpp>
#include <iostream>
using namespace cv;
using namespace std;
//低通滤波的类型：
enum LPFILTER_TYPE { ILP_FILTER = 0, BLP_FILTER = 1, GLP_FILTER = 2 };
namespace mycv
{
    void Im_Contrast(cv::Mat src,cv::Mat &dst);
    cv::Mat createLPFilter(Size size, Point center, float radius, int type, int n);
    //1-理想低通滤波器，2-巴特沃斯低通滤波器，3-高斯低通滤波器
    void fftFilter(cv::Mat srcImage,cv::Mat&dstImage,int rudus,int type);
    
}
//构造低通滤波器
cv::Mat mycv::createLPFilter(Size size, Point center, float radius, int type, int n)
{
    Mat lpFilter = Mat::zeros(size, CV_32FC1);
    int rows = size.height;
    int cols = size.width;
    if (radius <= 0)
        return lpFilter;
    //构造理想低通滤波器
    if (type == ILP_FILTER)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                float norm2 = pow(abs(float(r - center.y)), 2) + pow(abs(float(c - center.x)), 2);
                if (sqrt(norm2) < radius)
                    lpFilter.at<float>(r, c) = 1;
                else
                    lpFilter.at<float>(r, c) = 0;
            }
        }
    }
    //构造巴特沃斯低通滤波器
    if (type == BLP_FILTER)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                lpFilter.at<float>(r, c) = float(1.0 / (1.0 + pow(sqrt(pow(r - center.y, 2.0) + pow(c - center.x, 2.0)) / radius, 2.0*n)));
            }
        }
    }
    //构造高斯低通滤波
    if (type == GLP_FILTER)
    {
        for (int r = 0; r < rows; r++)
        {
            for (int c = 0; c < cols; c++)
            {
                lpFilter.at<float>(r, c) = float(exp(-(pow(c - center.x, 2.0) + pow(r - center.y, 2.0)) / (2 * pow(radius, 2.0))));
            }
        }
    }
    return lpFilter;
}
void mycv::Im_Contrast(cv::Mat src, cv::Mat & dst)
{
    int r = src.rows;
    int c = src.cols;
    double minValue;
    double maxValue;
    cv::Point minLocation;
    cv::Point maxLocation;
    cv::Point matchLocation;
    minMaxLoc(src, &minValue, &maxValue, &minLocation, &maxLocation, cv::Mat());
    src.convertTo(dst, CV_32FC1);
    cv::Mat src_32fc = dst + 0.001;
    cv::log(src_32fc, src_32fc);
    float src_aver =exp(cv::sum(src_32fc)[0] / (r*c));//平均值
    dst = (dst / src_aver + 1);
    cv::log(dst, dst);
    float m_denormitor;
    m_denormitor=cv::log(maxValue / src_aver + 1);
    dst = dst / m_denormitor;
    dst.convertTo(dst, CV_8UC1,255,0);
}

void mycv::fftFilter(cv::Mat srcImage, cv::Mat & dstImage, int radius, int type)
{
    //数据类型转换，转换为 浮点型
     
    Mat srcFc1;
    srcImage.convertTo(srcFc1, CV_32FC1, 1.0, 0.0);
    /* -- 第二步：每一个数乘以(-1)^(r+c) -- */
    for (int r = 0; r < srcFc1.rows; r++)
    {
        for (int c = 0; c < srcImage.cols; c++)
        {
            if ((r + c) % 2)
                srcFc1.at<float>(r, c) *= -1;
        }
    }
    //【2】将输入图像延扩到最佳的尺寸，边界用0补充
    int m = getOptimalDFTSize(srcImage.rows);
    int n = getOptimalDFTSize(srcImage.cols);
    //将添加的像素初始化为0.
    Mat padded;
    copyMakeBorder(srcFc1, padded, 0, m - srcImage.rows, 0, n - srcImage.cols, BORDER_CONSTANT, Scalar::all(0));
    //【3】为傅立叶变换的结果(实部和虚部)分配存储空间。
    //将planes数组组合合并成一个多通道的数组complexI
    Mat planes[] = { Mat_<float>(padded), Mat::zeros(padded.size(), CV_32F) };
    Mat complexI;
    merge(planes, 2, complexI);
    //【4】进行就地离散傅里叶变换
    dft(complexI, complexI);
    //【5】将复数转换为幅值，即=> log(1 + sqrt(Re(DFT(I))^2 + Im(DFT(I))^2))
    split(complexI, planes); // 将多通道数组complexI分离成几个单通道数组，planes[0] = Re(DFT(I), planes[1] = Im(DFTI))
    magnitude(planes[0], planes[1], planes[0]);// planes[0] = magnitude  
    Mat magnitudeImage = planes[0];//幅度谱
    //【6】进行对数尺度(logarithmic scale)缩放
    magnitudeImage += Scalar::all(1);
    log(magnitudeImage, magnitudeImage);//求自然对数
    normalize(magnitudeImage, magnitudeImage, 0, 1, NORM_MINMAX);
    //为了进行灰度级显示，做类型转换
    //magnitudeImage.convertTo(magnitudeImage, CV_8UC1, 255, 0);
    //找到傅里叶谱的最大值的坐标
    cv::Point maxLoc;//傅里叶谱的最大值的坐标
    minMaxLoc(magnitudeImage, NULL, NULL, NULL, &maxLoc);

    cv::Mat lpFilter = createLPFilter(magnitudeImage.size(), maxLoc, radius, type, 2);

    cv::Mat F_lpFilter;

    F_lpFilter.create(complexI.size(), complexI.type());

    for (int r = 0; r < F_lpFilter.rows; r++)
    {
        for (int c = 0; c < F_lpFilter.cols; c++)
        {
            //分别取出当前位置的快速傅里叶变换和理想低通滤波器的值
            Vec2f F_rc = complexI.at<Vec2f>(r, c);
            float lpFilter_rc = lpFilter.at<float>(r, c);
            //低通滤波器和图像的快速傅里叶变换对应位置相乘
            F_lpFilter.at<Vec2f>(r, c) = F_rc*lpFilter_rc;
        }
    }

    cv::dft(F_lpFilter, dstImage, DFT_SCALE + DFT_INVERSE + DFT_REAL_OUTPUT);
    /* -- 第九步：同乘以(-1)^(x+y) -- */
    for (int r = 0; r < dstImage.rows; r++)
    {
        for (int c = 0; c < dstImage.cols; c++)
        {
            if ((r + c) % 2)
                dstImage.at<float>(r, c) *= -1;
        }
    }
    //注意将结果转换 CV_8U 类型
    normalize(dstImage, dstImage, 0,255, NORM_MINMAX);
    dstImage.convertTo(dstImage, CV_8UC1);
    dstImage = dstImage(Rect(0, 0, srcImage.cols, srcImage.rows)).clone();

}