傅立葉轉換(Fourier Transform)

傅立葉轉換是一對一函數,可以是連續函數或者離散數列,正向傅立葉轉換,是把一個複雜的波,拆解成N個sin和cos組成的波,頻率從0倍到N-1倍,逆向傅立葉轉換,是把N個sin和cos組成的波,頻率從0倍到N-1倍,分別乘上強度、加上相位,再疊加成一個複雜的波。基本上任何的函式可以被無窮的sin和cos函式的加權和來表示,在影像處理上,經由傅立葉轉換將影像從空間域轉成頻率域,經過一些處理,再由反傅立葉轉換從頻率域轉回空間域。

傅立葉轉換時間複雜度為O(N2),實務上通常使用快速傅立葉轉換(Fast Fourier Transform, FFT),將公式的偶數項與奇數項分開整理,把時間複雜度降至O(NlogN),因為必須剛好對半分,所以影像的長、寬都須為2的次方,當長或寬不是2的次方,可在輸入像素末端補零,使得長和寬皆為2的次方。


Fourier Transform

以下示範將影像進行傅立葉轉換,得到頻譜,再從頻譜進行逆向傅立葉轉換,得到原始圖:

int main(){
    Mat inputImg = imread("lena.jpg", CV_LOAD_IMAGE_GRAYSCALE);
    Mat padded;                         
    int m = getOptimalDFTSize(inputImg.rows);  //m為大於等於inputImg.rows裡的最小值,且須為2、3、5的次方相乘
    int n = getOptimalDFTSize(inputImg.cols); 
    copyMakeBorder(inputImg, padded, 0, m-inputImg.rows, 0, n-inputImg.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);            

    split(complexImg, planes);                  //分離通道,planes[0]為實數部分,planes[1]為虛數部分 
    magnitude(planes[0], planes[1], planes[0]); //planes[0] = sqrt((planes[0])^2 + (planes[1])^2
    Mat magI = planes[0];
    magI += Scalar::all(1);                     //magI = log(1+planes[0])
    log(magI, magI);

    magI = magI(Rect(0, 0, magI.cols & -2, magI.rows & -2));  //令邊長為偶數

    //將區塊重排,讓原點在影像的中央
    int cx = magI.cols/2;
    int cy = magI.rows/2;

    Mat q0(magI, Rect(0, 0, cx, cy));   
    Mat q1(magI, Rect(cx, 0, cx, cy));  
    Mat q2(magI, Rect(0, cy, cx, cy));  
    Mat q3(magI, Rect(cx, cy, cx, cy)); 

    Mat tmp;                          
    q0.copyTo(tmp);
    q3.copyTo(q0);
    tmp.copyTo(q3);
    q1.copyTo(tmp);                    
    q2.copyTo(q1);
    tmp.copyTo(q2);

    normalize(magI, magI, 0, 1, CV_MINMAX); 

    imshow("輸入圖", inputImg);    
    imshow("頻譜", magI);

    //逆向傅立葉轉換
    Mat ifft;  
    idft(complexImg,ifft,DFT_REAL_OUTPUT);  
    normalize(ifft,ifft,0,1,CV_MINMAX);  
    imshow("逆向求輸入圖",ifft);  
    waitKey();

    return 0;    
}

Fourier Transform

Fourier Transform

Fourier Transform

繼續閱讀 傅立葉轉換(Fourier Transform)

SIFT特徵

在不同影像上進行特徵匹配時,常會遇到尺度變化的問題,也就是要分析的物體,可能在不同張影像的大小是不同的,當我們實際上要進行匹配時,由於尺度的差異,同個物體的特徵並不會匹配。

為了解決這個問題,有些算法用來尋找尺度不變的特徵,主要是基於每個檢測到的特徵點都伴隨著對應的尺寸,這邊介紹一種尺度不變的特徵,SIFT(Scale-Invariant Feature Transform)特徵,和另一個著名的尺度不變特徵檢測器SURF(Speeded Up Robust Features)相比,SURF速度較快,SIFT搜尋的特徵較準確。


以下使用程式碼使用SiftFeatureDetector找到影像的SURF特徵,接著用drawKeyPoints()劃出特徵點,旗標使用DRAW_RICH_KEYPOINTS,可看出圓圈的尺寸與特徵的尺度成正比,同時由於SURF算法將方向與每個特徵作關聯,使得特徵具有旋轉無關性:


影像匹配除了需要檢測圖像的特徵點,還需要向量來描述特徵才能進行比較,OpenCV裡使用SiftDescriptorExtractor來得到特徵點的SIFT描述子,描述子是一個Mat,行數與特徵點個數相同,每行都是一個N維的特徵向量,SIFT的描述子維度為128,特徵向量描述特徵點周圍的強度樣式,兩個特徵點越類似,特徵向量的距離越近。

假使我們想匹配在不同圖像的兩個相同物體,首先取得每個圖像的特徵,然後提取他們的描述子,將第一幅圖像和第二幅圖像中的每個特徵向量做比較,距離最近的特徵向量即為那個特徵的最佳匹配,對第一幅圖像的所有特徵都進行此處理,這也是BruteForceMatcher所採用的方法。

以下範例使用SIFT演算法檢測特徵點,並呼叫drawMatches()看特徵點的匹配情形:

SURF特徵

在不同影像上進行特徵匹配時,常會遇到尺度變化的問題,也就是要分析的物體,可能在不同張影像的大小是不同的,當我們實際上要進行匹配時,由於尺度的差異,同個物體的特徵並不會匹配。

為了解決這個問題,有些算法用來尋找尺度不變的特徵,主要是基於每個檢測到的特徵點都伴隨著對應的尺寸,這邊介紹一種尺度不變的特徵,SURF(Speeded Up Robust Features)特徵,不僅具有尺度和旋轉的不變性,而且非常高效,和另一個著名的尺度不變特徵檢測器SIFT(Scale-Invariant Feature Transform))相比,SURF速度較快,SIFT搜尋的特徵較準確。


以下使用程式碼使用SurfFeatureDetector找到影像的SURF特徵,接著用drawKeyPoints()劃出特徵點,旗標使用DRAW_RICH_KEYPOINTS,可看出圓圈的尺寸與特徵的尺度成正比,同時由於SURF算法將方向與每個特徵作關聯,使得特徵具有旋轉無關性:


影像匹配除了需要檢測圖像的特徵點,還需要向量來描述特徵才能進行比較,OpenCV裡使用SurfDescriptorExtractor來得到特徵點的SURF描述子,描述子是一個Mat,行數與特徵點個數相同,每行都是一個N維的特徵向量,SURF的默認維度為64,特徵向量描述特徵點周圍的強度樣式,兩個特徵點越類似,特徵向量的距離越近。

假使我們想匹配在不同圖像的兩個相同物體,首先取得每個圖像的特徵,然後提取他們的描述子,將第一幅圖像和第二幅圖像中的每個特徵向量做比較,距離最近的特徵向量即為那個特徵的最佳匹配,對第一幅圖像的所有特徵都進行此處理,這也是BruteForceMatcher所採用的方法。

以下範例使用SURF演算法檢測特徵點,並呼叫drawMatches()看特徵點的匹配情形:

FAST特徵

Harris算法提出了搜尋角點(或者說是特徵點)的方式,基於兩個正交方向上的強度變化率,能夠得到不錯的結果,但由於需要耗時的計算影像的一階微分,而特徵點檢測通常只是整個複雜算法的第一步,所以這邊另外介紹FAST檢測器,正如名字所示,此算法能快速進行檢測,依賴較少像素來確定是否為特徵點。

FAST(Features from Accelerated Segment Test)算法,和Harris算法同樣要定義甚麼是角點,FAST的角點定義基於和周圍像素的強度關係,檢查候選像素周圍一圈的像素,與中心點差異較大的像素如果組成連續的圓弧,並且弧長大於圓周長的3/4,那麼就把這個點視為特徵點。

FAST算法使用額外的技巧進行加速,一開始先測試周圍圓上被90分割的四個點,比如說測試點的上、下、左、右四個像素,因為滿足FAST對角點的定義,則這四個點至少有三個和中心點的強度差異須大於閾值,實務上大部分的像素都可用這方法進行移除,因此使用上非常高效,原則上測試上的半徑應該是一個參數,實際上半徑為3可以兼顧結果及效率。

和Harris找角點相同,可以在找到的角點上執行非極大值抑制,而由於FAST算法能快速的檢測特徵點,在對執行速度要求的環境下,比如在高FPS的影片上進行視覺跟蹤時,此時可以考慮用FAST來搜尋特徵點。

Harris 角點

在影像中檢測特徵點時,角點可以做為一個重要的參考,因為角點是兩條邊緣的交點處,可以被精確定位,這和位於相同強度的區域不同,與物體輪廓的點也不同,輪廓點難以在其他影像的相同物體進行精確定位。

Harris特徵檢測器是一個經典的角點檢測方法,OpenCV使用cornerHarris()實現Harris角點偵測演算法,輸出結果為浮點數類型的影像,每個像素值為相對位置的角點強度,之後再用閾值進行二值化即可得到角點。


OpenCV Harris角點檢測

void cornerHarris(InputArray src, OutputArray dst, int blockSize, int ksize, double k, int borderType=BORDER_DEFAULT)

  • src:輸入圖,8位元或浮點數單通道圖。
  • dst:輸出圖,儲存Harris檢測結果,型態為CV_32FC1,尺寸和輸入圖相同。
  • blockSize:相鄰像素的尺寸。
  • ksize:Sobel算子的濾波器模板大小。
  • k:Harris參數,即為下面方程式的k值。
  • borderType:邊緣擴充方式。

為了偵測影像中的角點,Harris觀察一個假定點周圍小窗口內強度差的平方和,窗口大小為cornerHarris()的第三個參數blockSize,因為無法確定高強度變化的方向,因此在所有可能的方向計算,過程首先獲得強度變化最大的方向,接著檢查它垂直方向的變化是否也很強烈,同時滿足的話便是一個角點。

小窗口強度差平方和:

Harris

用泰勒展開式對算式進行近似:

Harris

對算式化簡:

eHarris

轉換成矩陣型式:

Harris

矩陣M是一個斜方差(Covariance)矩陣,代表所有方向上的強度變化率,矩陣內的一階微分通常是Sobel算子計算結果,cornerHarris()的第四個參數ksize為Sobel的模板尺寸,斜方差矩陣的特徵值,代表最大強度變化以及和它垂直的方向,如果這兩個特徵值都低,代表此點位於變化不大的區域,要是其中一個較高另一個較低,代表此點位於邊上,如果兩個特徵值都較高,代表此點位於角點上。所以我們尋找角點的方式就是擁有超過閾值的斜方差矩陣特徵值:

Harris Harris

用下式驗證兩個特徵值是否足夠高,當兩個特徵值都高時,此計算結果也高,這也是cornerHarris()在每個位置得到的分數,k的值為cornerHarris()的第五個參數,實務上通常在0.05~0.5之間能得到滿意的結果: Harris

以下我們用cornerHarris()得到harris的分數,接著自定義閾值求出原始圖的角點:

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

int main(){
    Mat src = imread("church.jpg",CV_LOAD_IMAGE_GRAYSCALE);
    Mat harrisStrength;
    cornerHarris(src, harrisStrength, 3, 3, 0.01);
    Mat corners;
    double thres = 0.0001;
    threshold(harrisStrength, corners, thres, 255, THRESH_BINARY);

    imshow("原始圖", src);
    imshow("角點圖", corners);
    waitKey(0);
    return 0;
}

由上面結果可看出,實際呼叫cornerHarris()找角點時,偵測的結果可能在鄰近區域內有許多角點,而不容易進行精確定位,以下程式碼對Harris結果進行改進,角點不只需要結果高於閾值,還必須為局部最大值。所以對Harris結果圖進行膨脹,膨脹運算只有在局部最大值的地方維持原值,之後輸出結果為維持原值的位置才是角點,以下為詳細的程式碼:


關於特徵點聚集的問題,除了用局部極大值的方式,也可以指定兩個特徵點的最小距離,OpenCV另外有goodFeaturesToTrack(),從Harris得分最高的點開始,僅接受距離大於最小允許距離的特徵點,檢測的結果可用於視覺跟蹤的特徵集合。

相機校正(Camera calibration)

這邊介紹如何使用OpenCV提供的函式,對照相機進行校正。校正時我們必須先準備圖表,並用照相機在不同視角和距離對這個圖表拍照,使用時大約拍攝10到20張,以這些影像當作我們的輸入資料,目前OpenCV的支持三種類型的校準圖表,分別是:1.經典的黑白色棋盤,2.對稱圓形圖案,3.非對稱圓形圖案。

OpenCV在作相機校正時,有考量到徑向(radial)和切向(tangential)兩個因素,徑向失真會產生桶狀或類似魚眼效果的失真,對於一個校正前的像素點(X,Y),校正後的位置在(Xcorrected,Ycorrected),以下為徑向失真校正前後位置的關係: Camera calibration

切向失真主要發生在鏡頭和拍攝物體面不是平行,由以下的公式進行切向失真校正: Camera calibration

綜合以上兩個公式,OpenCV主要考量以下五個參數的矩陣進行校正: Camera calibration


OpenCV 找棋盤角點

bool findChessboardCorners(InputArray image, Size patternSize, OutputArray corners, int flags = CALIB_CB_ADAPTIVE_THRESH+CALIB_CB_NORMALIZE_IMAGE)

  • image:輸入圖,必須為8位元的灰階或彩色影像。
  • patternSize:棋盤的尺寸,patternSize = Size(points_per_row,points_per_colum) = Size(columns,rows)。
  • corners:輸出角點。
  • flags:旗標,CV_CALIB_CB_ADAPTIVE_THRESH表示使用區域自適應濾波進行二值化,CV_CALIB_CB_NORMALIZE_IMAGE表示二值化前先呼叫equalizeHist()。
  • 有找到棋盤角點返回true,否則返回false。

OpenCV 找精確角點

void cornerSubPix(InputArray image, InputOutputArray corners, Size winSize, Size zeroZone, TermCriteria criteria)

  • image:輸入圖。
  • corners:輸入角點,會更新成更精確的位置。
  • winSize:搜尋範圍,假設輸入Size(5,5),則會搜尋5*2+1的範圍,也就是11×11的尺寸
  • criteria:迭代搜尋的終止條件。

OpenCV 相機校正

double calibrateCamera(InputArrayOfArrays objectPoints, InputArrayOfArrays imagePoints, Size imageSize, InputOutputArray cameraMatrix, InputOutputArray distCoeffs, OutputArrayOfArrays rvecs, OutputArrayOfArrays tvecs, int flags=0, TermCriteria criteria=TermCriteria( TermCriteria::COUNT+TermCriteria::EPS, 30, DBL_EPSILON))

  • objectPoints:校正座標系統的位置,輸入型態為std::vector< std::vector< cv::Vec3f > >。
  • imagePoints:校正點的投影,輸入型態為std::vector< std::vector< cv::Vec2f > >,imagePoints.size()和objectPoints.size()必須相同,且imagePoints[i].size()和objectPoints[i].size()也必須相同。
  • imageSize:影像尺寸。
  • cameraMatrix:輸出3×3浮點數的相機矩陣。
  • distCoeffs:輸出的畸變參數,(k_1, k_2, p_1, p_2[, k_3[, k_4, k_5, k_6]]) 的4、5或8個元素。

OpenCV 校正矩陣

void initUndistortRectifyMap(InputArray cameraMatrix, InputArray distCoeffs, InputArray R, InputArray newCameraMatrix, Size size, int m1type, OutputArray map1, OutputArray map2)

  • cameraMatrix:輸入的相機矩陣。
  • distCoeffs:輸入的畸變參數。
  • newCameraMatrix:新的相機矩陣。
  • size:沒有畸變影像的尺寸。
  • m1type:map1的型態,可以為CV_32FC1或CV_16SC2。
  • map1:第一個輸出矩陣。
  • map2:第二個輸出矩陣。

我們程式碼用以下步驟進行相機校正:

  1. 多次拍攝圖表,並將影像存在硬碟作為輸入資料,這邊選用14張圖表影像。
  2. 尋找每張影像的校正點,為了得到更精確的位置,另外呼叫cornerSubPix()函式,TermCriteria 代表終止條件,這邊是用迭代次數30次,最小精度0.1當作終止條件。
  3. 輸入每張影像校正後的校正點位置。
  4. 假設這張影像的每個校正點都有找到,依我們的輸入圖來說,也就是24個點都有找到的話,就把這張影像的校正點存入當作輸入的資料點。
  5. 呼叫calibrateCamera()得到相機矩陣cameraMatrix,以及畸變矩陣distCoeffs,這兩個矩陣會在後續校正時使用。
  6. 呼叫initUndistortRectifyMap()得到映射位置矩陣map1、map2。
  7. 得到map1和map2後,即可對實際輸入的影像進行映射得到校正後的影像。
  8. 可到此處下載圖檔:下載圖檔

calibration.h

#include <cstdio>
#include <opencv2/opencv.hpp>
#include <vector>
#include <string>

using namespace cv;
using namespace std;

class CameraCalibrator{
private:
    vector<string> m_filenames;
    Size m_borderSize;
    vector<vector<Point2f> > m_srcPoints;
    vector<vector<Point3f> > m_dstPoints;
public:
    void setFilename();
    void setBorderSize(const Size &borderSize);
    void addChessboardPoints();
    void addPoints(const vector<Point2f> &imageCorners, const vector<Point3f> &objectCorners);
    void calibrate(const Mat &src, Mat &dst);
};

calibration.cpp

#include "calibration.h"

void CameraCalibrator::setFilename(){
    m_filenames.clear();
    m_filenames.push_back("chessboard01.jpg"); 
    m_filenames.push_back("chessboard02.jpg"); 
    m_filenames.push_back("chessboard03.jpg"); 
    m_filenames.push_back("chessboard04.jpg"); 
    m_filenames.push_back("chessboard05.jpg"); 
    m_filenames.push_back("chessboard06.jpg"); 
    m_filenames.push_back("chessboard07.jpg"); 
    m_filenames.push_back("chessboard08.jpg"); 
    m_filenames.push_back("chessboard09.jpg"); 
    m_filenames.push_back("chessboard10.jpg"); 
    m_filenames.push_back("chessboard11.jpg"); 
    m_filenames.push_back("chessboard12.jpg"); 
    m_filenames.push_back("chessboard13.jpg"); 
    m_filenames.push_back("chessboard14.jpg");
} 

void CameraCalibrator::setBorderSize(const Size &borderSize){ 
    m_borderSize = borderSize; 
} 

void CameraCalibrator::addChessboardPoints(){  
    vector<Point2f> srcCandidateCorners; 
    vector<Point3f> dstCandidateCorners; 
    for(int i=0; i<m_borderSize.height; i++){ 
        for(int j=0; j<m_borderSize.width; j++){ 
            dstCandidateCorners.push_back(Point3f(i, j, 0.0f)); 
        } 
    } 

    for(int i=0; i<m_filenames.size(); i++){ 
        Mat image=imread(m_filenames[i],CV_LOAD_IMAGE_GRAYSCALE); 
        findChessboardCorners(image, m_borderSize, srcCandidateCorners); 
        TermCriteria param(TermCriteria::MAX_ITER + TermCriteria::EPS, 30, 0.1);
        cornerSubPix(image, srcCandidateCorners, Size(5,5), Size(-1,-1), param);  
        if(srcCandidateCorners.size() == m_borderSize.area()){ 
            addPoints(srcCandidateCorners, dstCandidateCorners); 
        } 
    } 
} 

void CameraCalibrator::addPoints(const vector<Point2f> &srcCorners, const vector<Point3f> &dstCorners){           
    m_srcPoints.push_back(srcCorners);  
    m_dstPoints.push_back(dstCorners); 
} 

void CameraCalibrator::calibrate(const Mat &src, Mat &dst){ 
    Size imageSize = src.size(); 
    Mat cameraMatrix, distCoeffs, map1, map2; 
    vector<Mat> rvecs, tvecs; 
    calibrateCamera(m_dstPoints, m_srcPoints, imageSize, cameraMatrix, distCoeffs, rvecs, tvecs);  
    initUndistortRectifyMap(cameraMatrix, distCoeffs, Mat(), Mat(), imageSize, CV_32F, map1, map2); 
    remap(src, dst, map1, map2, INTER_LINEAR); 
}

main.cpp

#include "calibration.h"

int main(){ 
    CameraCalibrator myCameraCalibrator; 
    myCameraCalibrator.setFilename(); 
    myCameraCalibrator.setBorderSize(Size(6,4)); 
    myCameraCalibrator.addChessboardPoints();

    Mat src = imread("chessboard09.jpg",0);
    Mat dst;
    myCameraCalibrator.calibrate(src, dst);

    imshow("Original Image", src);
    imshow("Undistorted Image", dst);
    waitKey();
    return 0;
}

Camera calibration

Camera calibration

繼續閱讀 相機校正(Camera calibration)

GPU平行運算

GPU(graphics processing unit)一開始是為渲染圖形場景而建立,這些場景的數據不需要以串列的方式一步一步執行,而是用並行的方式一次性渲染大量的數據,從結構上來說,GPU不像CPU基於數個寄存器和高速指令集,GPU一般有數百個較小的處理單元。這些處理單元每一個都比CPU的核心慢很多,然而由於它的數目眾多,能夠同時進行大量運算,已經有越來越多的人嘗試用GPU來進行圖形渲染以外的工作,這就產生了GPGPU(general-purpose computation on graphics processing units)的概念。

圖像處理器有它自己的內存,一般稱呼為顯存,當我們從硬碟讀數據並產生一個Mat對象的時候,數據是放在普通內存當中的(由CPU掌管),CPU可以直接操作內存,然而GPU不能這樣,它只能操作它自己的顯存,所以需要先讓CPU將用於計算的信息轉移到GPU掌管的顯存上。完成這一個上傳過程,需要比訪問內存多很多的時間,而且最終的計算結果需要回傳到系統內存,由於傳輸的代價高昂,所以耗時太少的函數到GPU上計算只會造成反效果。

Mat對象僅僅存儲在內存中,GPU上的計算必須使用GpuMat,它的工作方式類似Mat,唯一的限制是不能直接引用GPU函數,使用時要傳輸Mat對像到​GPU上,且計算時使用gpu::內的相對函式,計算完成需再回傳結果,可使用簡單的重載操作符號或者調用下載函數。


以下我們比較傳統的