圖像的二維頻譜圖的理解 20170622

本文轉載自查看原文 2017-07-29 15:24 1262 零碎的數學知識點/ Matlab

# 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++ ->庫目錄）路徑：MATLAB2013a\extern\lib\win64\microsoft
// 
// 作者：路邊的十元錢硬幣
// 時間：20170630
// 最后的話，做注釋可以的話，最好還是英文，因為vs平台，漢字注釋有可能會造成編譯錯誤。
//

#include "stdafx.h"
#include <opencv2\opencv.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)的圖像）