• 图像的二维频谱图的理解 20170622


    # 1 图像二维频谱长什么样子(左图是原图,右图是对应的频谱图)

     

    (图片来源:第一组是来自matlab自带的图片 “cameraman.tif”;第二组是用 excel 画的,然后截图)

    # 2 怎么获得(matlab和C++调用)

    • matlaba代码,保存为 spectrum2D.m 
    function [Result] = spectrum2D(I)
    % I is a gray image
    % 'Input Image should be gray!' A = rgb2gray(I);
    
    A = I; % A 要是一个灰度图像
    F = fftshift(fft2(A)); % fft2() 是快速傅里叶变换函数
                           % fftshift() 目的是为了将图像横纵划分成四块,对角线互换。
                           % 经过变换以后的图像 F 中心是低频,往外频率增大。
    
    % 以下的操作是为获得可以显示的图像。
    % 如果直接显示 F ,得到的很可能是一个全白的图。
    % 原因是,F内像素的值一般都很大,远大于255.
    % 而matlab的imshow函数黑白对应0-255,所以不做归一化处理的图像,显示都是白色的。
    
    B = abs(real(F)); % F 有可能是负值,而且数据是中心对称的,中心为最亮点(可以参考C++例程中的说明)
    [m, n] = size(B); % 获得行数、列数
    R = reshape(B, 1, m*n); % 将B映射成1行m*n列的矩阵
    Temp = mapminmax(R, 0, 255); % 将R的像素值归一化到0-255之间
    Result = reshape(Temp, m, n); 
    
    % imshow(Result); % 建议调用的时候使用:imshow(spectrum2D(灰度图像矩阵))
                      % imshow(spectrum2D( rgb2gray(RGB图像矩阵)))
    
    end

    (以上三幅图,分别对应代码中的:A、F、Result)

    抱着折腾的心态,我们把matlab的这个函数在opencv里面试着来玩一下~~不感兴趣的可以跳过。

    使用opencv来调用matlab的目的在于可以方便的移植,而且方便我们可以数据进行处理;

      比如这后面要做的峰点滤波,使得亮点更加明显。(为什么要这样做呢,那就要理解亮点的含义才可以,接着往下看。)

    主要思路:1 使用matlab生成dll;2 在opencv环境中调用dll

    • matlab生成dll
    mex -setup
    mbuild -setup
    # mcc -W cpplib:名字  -T link:lib matlab文件.m
    mcc -W cpplib:spectrum2D -T link:lib spectrum2D.m
    • 调用opencv
    // Spectrum2D.cpp : 定义控制台应用程序的入口点。
    // 配置手顺:
    // 1. vs2015(工程属性设置为64bit)+opencv320,64 bit windows,配置好(不懂得话去网上搜吧,很容易)。
    // 2. 加载matlab生产的dll、头文件、lib,不需要cpp。(不懂得话,百度搜“怎么加载dll”)
    // 3. 上步的头文件应该会提示 mclmcrrt.h mclcppclass.h找不到,这是matlab安装目录下的文件,加载到(属性->VC++ ->包含目录)
    // 4. 加载matlab库文件 mclmcrrt.lib mclmcr.lib
    //      (属性->VC++ ->库目录)路径:MATLAB2013aexternlibwin64microsoft
    // 
    // 作者:路边的十元钱硬币
    // 时间:20170630
    // 最后的话,做注释可以的话,最好还是英文,因为vs平台,汉字注释有可能会造成编译错误。
    //
    
    #include "stdafx.h"
    #include <opencv2opencv.hpp>
    #include "spectrum2D.h"
    
    /// mwArray 是matlab用的数据类型,不论是矩阵还是数字,当需要作为参数传递到matlab生产的函数中时,要转化成这种数据类型
    /// 使用方法见代码
    mwArray Mat2mwArray(cv::Mat src)
    {
        CV_Assert(src.type() == CV_8UC1);
    
        mwArray dst(src.rows, src.cols, mxUINT8_CLASS); /// 初始化,可以理解成矩阵
        cv::Mat src_t = src.t();
        dst.SetData(src_t.data, src.rows*src.cols); /// 赋值
    
        return dst;
    }
    
    cv::Mat mwArry2Mat(mwArray src, int rows, int cols)
    {
        if(src.IsEmpty()) /// 是否为空
            return cv::Mat();
    
        cv::Mat dst = cv::Mat::zeros(rows, cols, CV_64FC1);
        for(int j(0); j<rows; ++j)
        {
            double* pdata = dst.ptr<double>(j);
            for(int i(0); i<cols; ++i)
            {
                pdata[i] = src(j+1,i+1); /// 元素访问(行号,列号)
            }
        }
        
        return dst;    
    }
    
    
    
    int _tmain(int argc, _TCHAR* argv[])
    {
        cv::Mat src = cv::imread("1.JPG", CV_LOAD_IMAGE_GRAYSCALE);
        cv::imshow("src", src);
    
        /// matlab mwArray 初始化,必须要做,不然报错
        if( !spectrum2DInitialize())
        {
            std::cout << "Could not initialize!" << std::endl;
            return -1; 
        }
    
        mwArray mlImg = Mat2mwArray(src);
        mwArray mlResult ;
    
        spectrum2D(1,mlResult, mlImg);
    
        /**
        简单说下mwArray的用法:
        mwArray a_int_varible(1, 1, mxINT8_CLASS); /// 定义
        mwArray a_double_varible(1, 1, mxDOUBLE_CLASS);
    
        a_int_varible(1,1) = 10; /// 赋值
        a_double_varible(1,1) = 10.0;
    
        int a = a_int_varible(1,1);
    
        **/
    
        cv::Mat result = mwArry2Mat(mlResult, src.rows, src.cols);
    
        cv::imshow("FFT2", result);
        cv::waitKey(0);
    
        return 0;
    }

    设置断点可以查看result矩阵中的数值。

    vs需要安装插件 Image Watch

     原图是#1中的第二组图,这是其对应频谱的中心亮点的局部放大图。其中255是中心亮点,6.8616是次亮点,注意左右是对称的。

    规律一:中心是最大亮点

    规律二:次亮点6.8616到255的距离是6,正好等于#1第二组图黑白条纹的重复次数。这不是巧合

    # 3 图的含义

    • 先看黑白条纹的图,把图像想象成波浪,将像素点值当成起伏高度,于是获得一个水平往右(往左也无所谓)传播的“黑色波浪”。
    • 在一维中,一条曲线可以被表示成很多不同频率的sin或者cos的叠加。(要问为什么的话,百度搜搜傅里叶变换相关的东西。)注意,这里的sin或者cos是二维波,只是线(y=sin(x)的图像)

    • 类似的,在图像中,每个“波浪”也可以被表示成很多不同频率的sin或者cos的叠加。注意,这里的sin或者cos函数是三维波,是曲面(y=sin(x),z=任意数)(波浪)

    • 想象一幅拥有最大频率波浪的图应该是什么样的,看图。

    是不是应该到最后,黑或者白条纹的宽度=1个像素?对的。这就是一幅图频率的极限。

    所以,得到重要结论:

    • 频谱图的中心亮点,是0频率点。往外,频率增大;同一圆周上的点,频率相同。
    • 频谱图的x轴的最右边点(无论是不是亮点),表示图像水平方向上的最大频率(波长是2个像素)。
    • 频谱图的任意点A到中心点O的距离|OA|表示频率

      大家应该是有疑惑的,就是这里的关于“频率”的理解。这里说到的“频率”,可以理解成条纹出现的次数,看下图说明。

        上图蓝色的为频谱图区域。其上的点A到中心点O的距离|AO|(单位像素),表示波的频率w=|OA|;其A点的像素值,表示图像中该方向(矢量OA的方向)上,频率为|AO|的比重(这里说比重而不说数量,因为频谱图的点的像素值只是响应值,只能用于对比多少。而且一般需要归一化。谈相差比例才有意义)

    • 总结:

    频谱图上任意点A,中心点O,A点像素值PA

    A点的含义是:原图像上,(矢量OA)方向上的,频率为|OA|的平面波的“比重”是PA。

    # 4 图的作用

    |OA| = 对应A点的波浪在原图上出现的次数。

    滤去某些噪声等。

    能力有限,错误或不全之处请不吝赐教。

    路边的十元钱硬币

    end.

  • 相关阅读:
    POJ Problem 1363 Rails 【栈】
    POJ Problem 3040 Allowance 【贪心】
    ACM 程序对拍
    HDU Problem
    POJ
    HDU Problem
    HDU Problem—2124 Repair the Wall 【贪心】
    HDU Problem 1052 Tian Ji -- The Horse Racing 【贪心】
    POJ Problem Radar Installation 【贪心】
    Beyond Compare和输出文件比较的方法
  • 原文地址:https://www.cnblogs.com/alexYuin/p/7067381.html
Copyright © 2020-2023  润新知