如何識別出輪廓准確的長和寬

問題來源：

實際項目中，需要給出識別輪廓的長度和寬度。

初步分析：

輪廓分析的例程為：

int main( int argc, char * * argv )
{
    //read the image
    Mat img = imread( "e:/sandbox/leaf.jpg");
    Mat bw;
    bool dRet;
    //resize
    pyrDown(img,img);
    pyrDown(img,img);

    cvtColor(img, bw, COLOR_BGR2GRAY);
    //morphology operation  
    threshold(bw, bw, 150, 255, CV_THRESH_BINARY);
    //bitwise_not(bw,bw);
    //find and draw contours
    vector <vector <Point > > contours;
    vector <Vec4i > hierarchy;
    findContours(bw, contours, hierarchy, CV_RETR_LIST, CV_CHAIN_APPROX_NONE);
    for ( int i = 0;i <contours.size();i ++)
    {
        RotatedRect minRect = minAreaRect( Mat(contours[i]) );
        Point2f rect_points[ 4];
        minRect.points( rect_points ); 
        for( int j = 0; j < 4; j ++ )
            line( img, rect_points[j], rect_points[(j + 1) % 4],Scalar( 255, 255, 0), 2);
    }
    imshow( "img",img);
    waitKey();
    return 0;
}

得到結果：

 

 

對於這樣 的輪廓分析，標明出來的1和2明顯是錯誤的。但是除了minAreaRect之外，已經沒有更解近一步的方法。

也嘗試首先對輪廓進行凸包處理，再查找外接矩形，效果同樣不好。

 

解題思路：

仍然要從現有的、穩定運行的代碼里面找方法。目前OpenCV函數 getOrientation 能夠通過PCA方法找到圖像/輪廓的方向

比如這樣：

 

 

在項目圖片上能夠得到這樣結果：

顯然是更符合實際情況的，當然，葉柄這里產生了干擾，但那是另一個問題。

獲得主方向后，下一步就是如何獲得准確的長和寬。PCA方法無法獲得長寬，也嘗試通過旋轉矩陣的方法直接獲得結果：

////以RotatedRect的方式返回結果

    //RotatedRect box;

    //box.center.x = pos.x;

    //box.center.y = pos.y;

    //box.size.width = flong;

    //box.size.height = fshort;

    //box.angle = (float)atan2( eigen_vecs[0].y, eigen_vecs[0].x)*180/3.1415926; //弧度轉角度

 

    ////繪制rotateRect

    //Point2f rect_points[4];

    //box.points( rect_points ); 

    //for( int j = 0; j < 4; j++ )

    //    line( img, rect_points[j], rect_points[(j+1)%4],Scalar(0,0,255),2);

但是需要注意的是，這里的pca獲得的center並不是絕對的center，而且在中線兩邊，輪廓到中線的長度不一定一樣。為了獲得最精確的結果，就需要直接去求出每個邊的長度，並且繪制出來。

思路很簡單，就是通過中線（及其中線的垂線）將原輪廓分為兩個部分，分別求這兩個部分的到中線的最大距離（加起來就是長，分開來就是位置）。

求的長軸端點：

 

 

求得到中線最遠距離點（藍色），這也就是到中線的距離。

 

 

 

 

距離的計算很多時候只是點的循環。最后存在一個問題，那就是這樣一個圖像，已經知道p0-03的坐標，和兩條軸線的斜率，如何繪制4個

角點？

 

實際上，這是一個數學問題，並且有解析解：

    //通過解析方法，獲得最后結果 

    Point p[4]; 

    p[0].x = (k_long * _p[0].x   - k_short * _p[2].x  +  _p[2].y - _p[0].y)  / (k_long - k_short);

    p[0].y = (p[0].x - _p[0].x)*k_long + _p[0].y;

    p[1].x = (k_long * _p[0].x   - k_short * _p[3].x  +  _p[3].y - _p[0].y)  / (k_long - k_short);

    p[1].y = (p[1].x - _p[0].x)*k_long + _p[0].y;

    p[2].x = (k_long * _p[1].x   - k_short * _p[2].x  +  _p[2].y - _p[1].y)  / (k_long - k_short);

    p[2].y = (p[2].x - _p[1].x)*k_long + _p[1].y;

    p[3].x = (k_long * _p[1].x   - k_short * _p[3].x  +  _p[3].y - _p[1].y)  / (k_long - k_short);

    p[3].y = (p[3].x - _p[1].x)*k_long + _p[1].y;

 

成功！！！

得到最后結果，這正是我想要得到的。但是由於算法穩定性方面和效率的考慮，還需要進一步增強。

  

p.s

重新翻了一下minarearect

cv : :RotatedRect cv : :minAreaRect( InputArray _points )
{
    CV_INSTRUMENT_REGION()

    Mat hull;
    Point2f out[ 3];
    RotatedRect box;

    convexHull(_points, hull, true, true);

    if( hull.depth() != CV_32F )
    {
        Mat temp;
        hull.convertTo(temp, CV_32F);
        hull = temp;
    }

    int n = hull.checkVector( 2);
    const Point2f * hpoints = hull.ptr <Point2f >();

    if( n > 2 )
    {
        rotatingCalipers( hpoints, n, CALIPERS_MINAREARECT, ( float *)out );
        box.center.x = out[ 0].x + (out[ 1].x + out[ 2].x) * 0. 5f;
        box.center.y = out[ 0].y + (out[ 1].y + out[ 2].y) * 0. 5f;
        box.size.width = ( float)std : :sqrt(( double)out[ 1].x *out[ 1].x + ( double)out[ 1].y *out[ 1].y);
        box.size.height = ( float)std : :sqrt(( double)out[ 2].x *out[ 2].x + ( double)out[ 2].y *out[ 2].y);
        box.angle = ( float)atan2( ( double)out[ 1].y, ( double)out[ 1].x );
    }
    else if( n == 2 )
    {
        box.center.x = (hpoints[ 0].x + hpoints[ 1].x) * 0. 5f;
        box.center.y = (hpoints[ 0].y + hpoints[ 1].y) * 0. 5f;
        double dx = hpoints[ 1].x - hpoints[ 0].x;
        double dy = hpoints[ 1].y - hpoints[ 0].y;
        box.size.width = ( float)std : :sqrt(dx *dx + dy *dy);
        box.size.height = 0;
        box.angle = ( float)atan2( dy, dx );
    }
    else
    {
        if( n == 1 )
            box.center = hpoints[ 0];
    }

    box.angle = ( float)(box.angle * 180 /CV_PI);
    return box;
}

 

那么，這個官方函數首先就把輪廓找了hull矩。這個當然對於很多問題都是好方法，很簡單直觀（我這里的方法就繁瑣很多）。但是忽視了一個重要問題：hull變換后，會丟失信息。顯然這就是結果不准確的原因。

 

我也在answeropencv上進行了咨詢，berak給出的comment是

maybe:

1. find principal axes of your shape (PCA)
2. rotate to upright (warp)
3. boundingRect() (image axis aligned)

這個顯然也不是很妥當，因為這個結果還需要rotate回去，這會很麻煩。

感謝閱讀至此，希望有所幫助。

 

2018年8月29日19:43:01

經過一段時間后反思這個項目，應該說這個算法有一定自創的元素在里面，但是由於這個問題比較小眾，所以即使是在answeropencv上面，參與討論的人也比較少，使得算法有很多不充分的地方在里面：最主要的問題就是在算法的后面部分，多次進行全輪廓循環，使得算法的效率降低。

那么，有沒有能夠提速的方法了？還是之前做的一些數學的研究和在answeropencv上 berak的comment 提醒了我，注意看：

 

這里，黑色的是原始的OpenCV的坐標系，紅色的是新求出來的坐標系，你花了那么大功夫去算交點，實際上，不如將這個圖像旋轉為正，將外界矩形算出來，然后再反方向旋轉回去。這樣的話思路就很清楚了，但是需要花一點修改的時間。

但是正是因為這里的思路比較清晰，所以代碼寫起來，比較流暢，很快我就得到了下面的結果：

 

#include "stdafx.h"

#include "opencv2/imgcodecs.hpp"

#include "opencv2/highgui.hpp"

#include "opencv2/imgproc.hpp"

#include "opencv2/photo.hpp"

 

 

using namespace std;

using namespace cv;

#define DEBUG FALSE

 

//獲得單個點經過旋轉后所在精確坐標

Point2f GetPointAfterRotate(Point2f inputpoint,Point2f center,double angle){

    Point2d preturn;

    preturn.x = (inputpoint.x - center.x)*cos(-angle) - (inputpoint.y - center.y)*sin(-angle)+center.x;

    preturn.y = (inputpoint.x - center.x)*sin(-angle) + (inputpoint.y - center.y)*cos(-angle)+center.y;

    return preturn;

}

Point GetPointAfterRotate(Point inputpoint,Point center,double angle){

    Point preturn;

    preturn.x = (inputpoint.x - center.x)*cos(-1*angle) - (inputpoint.y - center.y)*sin(-1*angle)+center.x;

    preturn.y = (inputpoint.x - center.x)*sin(-1*angle) + (inputpoint.y - center.y)*cos(-1*angle)+center.y;

    return preturn;

}

 

double getOrientation(vector<Point> &pts, Point2f& pos,Mat& img)

{

    //Construct a buffer used by the pca analysis

    Mat data_pts = Mat(pts.size(), 2, CV_64FC1);

    for (int i = 0; i < data_pts.rows; ++i)

    {

        data_pts.at<double>(i, 0) = pts[i].x;

        data_pts.at<double>(i, 1) = pts[i].y;

    }

 

    //Perform PCA analysis

    PCA pca_analysis(data_pts, Mat(), CV_PCA_DATA_AS_ROW);

 

    //Store the position of the object

    pos = Point2f(pca_analysis.mean.at<double>(0, 0),

        pca_analysis.mean.at<double>(0, 1));

 

    //Store the eigenvalues and eigenvectors

    vector<Point2d> eigen_vecs(2);

    vector<double> eigen_val(2);

    for (int i = 0; i < 2; ++i)

    {

        eigen_vecs[i] = Point2d(pca_analysis.eigenvectors.at<double>(i, 0),

            pca_analysis.eigenvectors.at<double>(i, 1));

 

        eigen_val[i] = pca_analysis.eigenvalues.at<double>(i,0);

    }

 

    // Draw the principal components

    //在輪廓/圖像中點繪制小圓

    //circle(img, pos, 3, CV_RGB(255, 0, 255), 2);

    ////計算出直線，在主要方向上繪制直線

    //line(img, pos, pos + 0.02 * Point2f(eigen_vecs[0].x * eigen_val[0], eigen_vecs[0].y * eigen_val[0]) , CV_RGB(255, 255, 0));

    //line(img, pos, pos + 0.02 * Point2f(eigen_vecs[1].x * eigen_val[1], eigen_vecs[1].y * eigen_val[1]) , CV_RGB(0, 255, 255));

    return atan2(eigen_vecs[0].y, eigen_vecs[0].x);

}

 

//程序主要部分

int main( int argc, char** argv )

{

    //讀入圖像，轉換為灰度

    Mat img = imread("e:/sandbox/leaf.jpg");

    pyrDown(img,img);

    pyrDown(img,img);

 

    Mat bw;

    bool dRet;

    cvtColor(img, bw, COLOR_BGR2GRAY);

    //閾值處理

    threshold(bw, bw, 150, 255, CV_THRESH_BINARY);

    //尋找輪廓

    vector<vector<Point> > contours;

    vector<Vec4i> hierarchy;

    findContours(bw, contours, hierarchy, CV_RETR_LIST, CV_CHAIN_APPROX_NONE);

 

    //輪廓分析，找到

    for (size_t i = 0; i < contours.size(); ++i)

    {

        //計算輪廓大小

        double area = contourArea(contours[i]);

        //去除過小或者過大的輪廓區域（科學計數法表示）

        if (area < 1e2 || 1e5 < area) continue;

        //繪制輪廓

        drawContours(img, contours, i, CV_RGB(255, 0, 0), 2, 8, hierarchy, 0);

        //獲得輪廓的角度

        Point2f* pos = new Point2f();

        double dOrient =  getOrientation(contours[i], *pos,img);

        //轉換輪廓,並獲得極值

        for (size_t j = 0;j<contours[i].size();j++)

            contours[i][j] = GetPointAfterRotate(contours[i][j],(Point)*pos,dOrient);

        Rect rect = boundingRect(contours[i]);//輪廓最小外接矩形

        RotatedRect rotateRect = RotatedRect((Point2f)rect.tl(),Point2f(rect.br().x,rect.tl().y),(Point2f)rect.br());

        //將角度轉換回去並繪圖

        Point2f rect_points[4];

        rotateRect.points( rect_points ); 

        for (size_t j = 0;j<4;j++)

            rect_points[j] = GetPointAfterRotate((Point)rect_points[j],(Point)*pos,-dOrient);

        for( size_t j = 0; j < 4; j++ )

            line( img, rect_points[j], rect_points[(j+1)%4],Scalar(255,255,0),2);

        //得出結果    

        char cbuf[255];

        double fshort = std::min(rect.width,rect.height);

        double flong  = std::max(rect.width,rect.height);

        sprintf_s(cbuf,"第%d個輪廓,長度%.2f,寬度%.2f像素\n",i,flong,fshort);

    }

    return 0;

}

 

 

這段代碼中值得一提的是

Point GetPointAfterRotate(Point inputpoint,Point center,double angle){

    Point preturn;

    preturn.x = (inputpoint.x - center.x)*cos(-1*angle) - (inputpoint.y - center.y)*sin(-1*angle)+center.x;

    preturn.y = (inputpoint.x - center.x)*sin(-1*angle) + (inputpoint.y - center.y)*cos(-1*angle)+center.y;

    return preturn;

}

這個函數是直接計算出某一個點在旋轉后位置，采用的是數學方法推到，應該算自己創的函數。很多時候，我們並不需要旋轉整個圖像，而只是要獲得圖像旋轉以后的位置。

 

反思小結：應該說當時answerOpenCV上就給出了正確的結果提示，但是由於那時我鑽在自己的算法里面，沒能夠接受新的想法；過去一段時間后回顧，才發現了更好的解決方法。

但是走彎路並不可怕，只有不斷、持續地思考，盡可能將現有的解決方法優化，才可能在面對新的問題的時候有更多的手段、更容易提出創造出“方便書寫、效果顯著”的算法。

此外，基礎能力非常重要，如果基礎不牢，在創建新算法 的時候會遇到更多的困難，畢竟：基礎不牢、地動山搖。

感謝閱讀至此、希望有所幫助。

 

來自為知筆記(Wiz)

附件列表