Shi-Tomasi運算元

就像Harrise運算元是在Moravec運算元的基礎上改進得到的一樣：
http://blog.csdn.net/chaipp0607/article/details/54692818
Shi-Tomasi運算元是在Harrise運算元的基礎上改進的，改進之處在於他們使用了不同的響應函式。

Harrise運算元的響應函式為：

Shi-Tomasi運算元的響應函式為：

opencv實現

opencv提供了goodFeaturesToTrack函式用來實現Shi-Tomasi運算元，其API函式介面為：

void cv::goodFeaturesToTrack( 
InputArray _image, 
OutputArray _corners,  
int
 
 maxCorners,
double qualityLevel, 
double minDistance,  
InputArray _mask, 
int blockSize,  
bool useHarrisDetector, double harrisK )  

第一個引數：8位或32位浮點型輸入影象，單通道
第二個引數：儲存檢測出的角點
第三個引數：角點數目最大值，如果實際檢測的角點超過此值，則只返回前maxCorners個強角點
第四個引數：角點的品質因子
第五個引數：對於初選出的角點而言，如果在其周圍minDistance範圍內存在其他更強角點，則將此角點刪除
第六個引數：指定感興趣區，如不需在整幅圖上尋找角點，則用此引數指定ROI
第七個引數：計算協方差矩陣時的視窗大小
第八個引數：指示是否使用Harris角點檢測，如不指定，則計算shi-tomasi角點
第九個引數：Harris角點檢測需要的k值

goodFeaturesToTrack原始碼位置在：
….\opencv\sources\modules\imgproc\src\featureselect.cpp

原始碼為：

void cv::goodFeaturesToTrack( InputArray _image, OutputArray _corners,  
                              int maxCorners, double qualityLevel, double minDistance,  
                              InputArray _mask, int
 
 blockSize,  
                              bool useHarrisDetector, double harrisK )  
{  
    //如果需要對_image全圖操作，則給_mask傳入cv::Mat()，否則傳入感興趣區域  
    Mat image = _image.getMat(), mask = _mask.getMat();    

    CV_Assert( qualityLevel > 0 && minDistance >= 0 && maxCorners >= 0 );  //對引數有一些基本要求  
    CV_Assert( mask.empty() || (mask.type() == CV_8UC1 && mask.size() == image.size()) );  

    Mat eig, tmp;   //eig儲存每個畫素協方差矩陣的最小特徵值，tmp用來儲存經膨脹後的eig  
    if( useHarrisDetector )  
        cornerHarris( image, eig, blockSize, 3, harrisK ); //blockSize是計算2*2協方差矩陣的視窗大小，sobel運算元視窗為3，harrisK是計算Harris角點時需要的值  
    else  
        cornerMinEigenVal( image, eig, blockSize, 3 );  //計算每個畫素對應的協方差矩陣的最小特徵值，儲存在eig中  

    double maxVal = 0;  
    minMaxLoc( eig, 0, &maxVal, 0, 0, mask );   //maxVal儲存了eig的最大值  
    threshold( eig, eig, maxVal*qualityLevel, 0, THRESH_TOZERO );  //閾值設定為maxVal乘以qualityLevel，大於此閾值的保持不變，小於此閾值的都設為0  

    //預設用3*3的核膨脹，膨脹之後，除了區域性最大值點和原來相同，其它非區域性最大值點被    
    //3*3鄰域內的最大值點取代，如不理解，可看一下灰度影象的膨脹原理    
    dilate( eig, tmp, Mat());  //tmp中儲存了膨脹之後的eig  

    Size imgsize = image.size();   

    vector<const float*> tmpCorners;  //存放粗選出的角點地址  

    // collect list of pointers to features - put them into temporary image   
    for( int y = 1; y < imgsize.height - 1; y++ )  
    {  
        const float* eig_data = (const float*)eig.ptr(y);  //獲得eig第y行的首地址  
        const float* tmp_data = (const float*)tmp.ptr(y);  //獲得tmp第y行的首地址  
        const uchar* mask_data = mask.data ? mask.ptr(y) : 0;  

        for( int x = 1; x < imgsize.width - 1; x++ )  
        {  
            float val = eig_data[x];  
            if( val != 0 && val == tmp_data[x] && (!mask_data || mask_data[x]) )  //val == tmp_data[x]說明這是區域性極大值  
                tmpCorners.push_back(eig_data + x);  //儲存其位置  
        }  
    }  

    //-----------此分割線以上是根據特徵值粗選出的角點，我們稱之為弱角點----------//  
    //-----------此分割線以下還要根據minDistance進一步篩選角點，仍然能存活下來的我們稱之為強角點----------//  

    sort( tmpCorners, greaterThanPtr<float>() );  //按特徵值降序排列，注意這一步很重要，後面的很多程式設計思路都是建立在這個降序排列的基礎上  
    vector<Point2f> corners;  
    size_t i, j, total = tmpCorners.size(), ncorners = 0;  

    //下面的程式有點稍微難理解，需要自己仔細想想  
    if(minDistance >= 1)    
    {  
         // Partition the image into larger grids  
        int w = image.cols;  
        int h = image.rows;  

        const int cell_size = cvRound(minDistance);   //向最近的整數取整  

    //這裡根據cell_size構建了一個矩形視窗grid(雖然下面的grid定義的是vector<vector>，而並不是我們這裡說的矩形視窗，但為了便於理解,還是將grid想象成一個grid_width * grid_height的矩形視窗比較好)，除以cell_size說明grid窗口裡相差一個畫素相當於_image裡相差minDistance個畫素，至於為什麼加上cell_size - 1後面會講  
        const int grid_width = (w + cell_size - 1) / cell_size;   
        const int grid_height = (h + cell_size - 1) / cell_size;  

        std::vector<std::vector<Point2f> > grid(grid_width*grid_height);  //vector裡面是vector，grid用來儲存獲得的強角點座標  

        minDistance *= minDistance;  //平方，方面後面計算，省的開根號  

        for( i = 0; i < total; i++ )     // 剛剛粗選的弱角點，都要到這裡來接收新一輪的考驗  
        {  
            int ofs = (int)((const uchar*)tmpCorners[i] - eig.data);  //tmpCorners中儲存了角點的地址，eig.data返回eig記憶體塊的首地址  
            int y = (int)(ofs / eig.step);   //角點在原影象中的行  
            int x = (int)((ofs - y*eig.step)/sizeof(float));  //在原影象中的列  

            bool good = true;  //先認為當前角點能接收考驗，即能被保留下來  

            int x_cell = x / cell_size;  //x_cell，y_cell是角點（y,x）在grid中的對應座標  
            int y_cell = y / cell_size;  

            int x1 = x_cell - 1;  // (y_cell，x_cell）的4鄰域畫素  
            int y1 = y_cell - 1;  //現在知道為什麼前面grid_width定義時要加上cell_size - 1了吧，這是為了使得（y,x）在grid中的4鄰域畫素都存在，也就是說(y_cell，x_cell）不會成為邊界畫素  
            int x2 = x_cell + 1;    
            int y2 = y_cell + 1;  

            // boundary check，再次確認x1,y1,x2或y2不會超出grid邊界  
            x1 = std::max(0, x1);  //比較0和x1的大小  
            y1 = std::max(0, y1);  
            x2 = std::min(grid_width-1, x2);  
            y2 = std::min(grid_height-1, y2);  

            //記住grid中相差一個畫素，相當於_image中相差了minDistance個畫素  
            for( int yy = y1; yy <= y2; yy++ )  // 行  
            {  
                for( int xx = x1; xx <= x2; xx++ )  //列  
                {  
                    vector <Point2f> &m = grid[yy*grid_width + xx];  //引用  

                    if( m.size() )  //如果(y_cell，x_cell)的4鄰域畫素，也就是(y,x)的minDistance鄰域畫素中已有被保留的強角點  
                    {                 
                        for(j = 0; j < m.size(); j++)   //當前角點周圍的強角點都拉出來跟當前角點比一比  
                        {  
                            float dx = x - m[j].x;  
                            float dy = y - m[j].y;  
               //注意如果(y,x)的minDistance鄰域畫素中已有被保留的強角點，則說明該強角點是在(y,x)之前就被測試過的，又因為tmpCorners中已按照特徵值降序排列（特徵值越大說明角點越好），這說明先測試的一定是更好的角點，也就是已儲存的強角點一定好於當前角點，所以這裡只要比較距離，如果距離滿足條件，可以立馬扔掉當前測試的角點  
                            if( dx*dx + dy*dy < minDistance )  
                            {                                                         
                good = false;  
                                goto break_out;  
                            }  
                        }  
                    }  
                }   // 列  
            }    //行  

            break_out:  

            if(good)  
            {  
                // printf("%d: %d %d -> %d %d, %d, %d -- %d %d %d %d, %d %d, c=%d\n",  
                //    i,x, y, x_cell, y_cell, (int)minDistance, cell_size,x1,y1,x2,y2, grid_width,grid_height,c);  
                grid[y_cell*grid_width + x_cell].push_back(Point2f((float)x, (float)y));  

                corners.push_back(Point2f((float)x, (float)y));  
                ++ncorners;  

                if( maxCorners > 0 && (int)ncorners == maxCorners )  //由於前面已按降序排列，當ncorners超過maxCorners的時候跳出迴圈直接忽略tmpCorners中剩下的角點，反正剩下的角點越來越弱  
                    break;  
            }  
        }  
    }  
    else    //除了畫素本身，沒有哪個鄰域畫素能與當前畫素滿足minDistance < 1,因此直接儲存粗選的角點  
    {  
        for( i = 0; i < total; i++ )  
        {  
            int ofs = (int)((const uchar*)tmpCorners[i] - eig.data);  
            int y = (int)(ofs / eig.step);   //粗選的角點在原影象中的行  
            int x = (int)((ofs - y*eig.step)/sizeof(float));  //在影象中的列  

            corners.push_back(Point2f((float)x, (float)y));  
            ++ncorners;  
            if( maxCorners > 0 && (int)ncorners == maxCorners )    
                break;  
        }  
    }  

    Mat(corners).convertTo(_corners, _corners.fixedType() ? _corners.type() : CV_32F);  

    /* 
    for( i = 0; i < total; i++ ) 
    { 
        int ofs = (int)((const uchar*)tmpCorners[i] - eig.data); 
        int y = (int)(ofs / eig.step); 
        int x = (int)((ofs - y*eig.step)/sizeof(float)); 

        if( minDistance > 0 ) 
        { 
            for( j = 0; j < ncorners; j++ ) 
            { 
                float dx = x - corners[j].x; 
                float dy = y - corners[j].y; 
                if( dx*dx + dy*dy < minDistance ) 
                    break; 
            } 
            if( j < ncorners ) 
                continue; 
        } 

        corners.push_back(Point2f((float)x, (float)y)); 
        ++ncorners; 
        if( maxCorners > 0 && (int)ncorners == maxCorners ) 
            break; 
    } 
*/  
}  

一個測試例子：

#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
using namespace cv;
using namespace std;
int main()
{
    cv::Mat srcImage = cv::imread("1.jpg");
    if (!srcImage.data)
    return -1;
    cv::Mat srcGray;
    cvtColor(srcImage, srcGray, CV_BGR2GRAY);
    vector<Point2f> vecCorners;
    double qualityLevel = 0.02;
    double minDistance = 13;
    int blockSize = 3;
    bool useHarrisDetector = false;
    double k = 0.05;
    int maxCorners = 40;
    int maxTrackbar = 80;
    cv::Mat resultImage = srcImage.clone();
    goodFeaturesToTrack(srcGray,
        vecCorners,maxCorners,
        qualityLevel,minDistance,
        Mat(),blockSize,useHarrisDetector,k);
    std::cout << "角點個數為:" << vecCorners.size() << endl;
    for (int i = 0; i < vecCorners.size(); i++)
    {
        circle(resultImage, vecCorners[i], 4, 
            Scalar(0, 255, 255), 2);
    }
    cv::imshow("srcImage", srcImage);
    cv::imshow("resultImage", resultImage);
    cv::imwrite("srcImage.jpg",srcImage);
    cv::imwrite("resultImage.jpg",resultImage);
    cv::waitKey(0);
    return(0);
}

原圖：

結果：