如何识别出轮廓准确的长和宽

轮廓分析的例程为：

1. int main( int argc, char** argv )

2. {

3.    //read the image

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

5.    Mat bw;

6.    bool dRet;

7.    //resize

8.    pyrDown(img,img);

9.    pyrDown(img,img);

11.    cvtColor(img, bw, COLOR_BGR2GRAY);

12.    //morphology operation  

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

14.    //bitwise_not(bw,bw);

15.    //find and draw contours

16.    vector<vector<Point> > contours;

17.    vector<Vec4i> hierarchy;

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

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

20.    {

21.        RotatedRect minRect = minAreaRect( Mat(contours[i]) );

22.        Point2f rect_points[4];

23.        minRect.points( rect_points );

24.        for( int j = 0; j < 4; j++ )

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

26.    }

27.    imshow("img",img);

28.    waitKey();

29.    return 0;

30. }

得到结果：

对于这样 的轮廓分析，标明出来的1和2明显是错误的。但是除了minAreaRect之外，已经没有更解近一步的方法。也尝试首先对轮廓进行凸包处理，再查找外接矩形，效果同样不好。

解题思路：仍然要从现有的、稳定运行的代码里面找方法。目前OpenCV函数getOrientation能够通过PCA方法找到图像/轮廓的方向。比如这样：

在项目图片上能够得到这样结果：

显然是更符合实际情况的，当然，叶柄这里产生了干扰，但那是另一个问题。获得主方向后，下一步就是如何获得准确的长和宽。PCA方法无法获得长宽，也尝试通过旋转矩阵的方法直接获得结果：

1. 以RotatedRect的方式返回结果

2.    //RotatedRect box;

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

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

5.    //box.size.width = flong;

6.    //box.size.height = fshort;

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

8.    绘制rotateRect

9.    //Point2f rect_points[4];

10.    //box.points( rect_points );

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

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

但是需要注意的是，这里的pca获得的center并不是绝对的center，而且在中线两边，轮廓到中线的长度不一定一样。为了获得最精确的结果，就需要直接去求出每个边的长度，并且绘制出来。思路很简单，就是通过中线（及其中线的垂线）将原轮廓分为两个部分，分别求这两个部分的到中线的最大距离（加起来就是长，分开来就是位置）。求的长轴端点：

求得到中线最远距离点（蓝色），这也就是到中线的距离。

距离的计算很多时候只是点的循环。最后存在一个问题，那就是这样一个图像，已经知道p0-03的坐标，和两条轴线的斜率，如何绘制4个角点？

实际上，这是一个数学问题，并且有解析解：

1. //通过解析方法，获得最后结果

2. Point p[4];

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

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

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

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

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

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

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

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

成功！！！

得到最后结果，这正是我想要得到的。但是由于算法稳定性方面和效率的考虑，还需要进一步增强

p.s 重新翻了一下minarearect

1. cv::RotatedRect cv::minAreaRect( InputArray _points )

2. {

3.    CV_INSTRUMENT_REGION()

5.    Mat hull;

6.    Point2f out[3];

7.    RotatedRect box;

9.    convexHull(_points, hull, true, true);

11.    if( hull.depth() != CV_32F )

12.    {

13.        Mat temp;

14.        hull.convertTo(temp, CV_32F);

15.        hull = temp;

16.    }

18.    int n = hull.checkVector(2);

19.    const Point2f* hpoints = hull.ptr<Point2f>();

21.    if( n > 2 )

22.    {

23.        rotatingCalipers( hpoints, n, CALIPERS_MINAREARECT, (float*)out );

24.        box.center.x = out[0].x + (out[1].x + out[2].x)*0.5f;

25.        box.center.y = out[0].y + (out[1].y + out[2].y)*0.5f;

26.        box.size.width = (float)std::sqrt((double)out[1].x*out[1].x + (double)out[1].y*out[1].y);

27.        box.size.height = (float)std::sqrt((double)out[2].x*out[2].x + (double)out[2].y*out[2].y);

28.        box.angle = (float)atan2( (double)out[1].y, (double)out[1].x );

29.    }

30.    else if( n == 2 )

31.    {

32.        box.center.x = (hpoints[0].x + hpoints[1].x)*0.5f;

33.        box.center.y = (hpoints[0].y + hpoints[1].y)*0.5f;

34.        double dx = hpoints[1].x - hpoints[0].x;

35.        double dy = hpoints[1].y - hpoints[0].y;

36.        box.size.width = (float)std::sqrt(dx*dx + dy*dy);

37.        box.size.height = 0;

38.        box.angle = (float)atan2( dy, dx );

39.    }

40.    else

41.    {

42.        if( n == 1 )

43.            box.center = hpoints[0];

44.    }

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

47.    return box;

48. }

那么，这个官方函数首先就把轮廓找了hull矩。这个当然对于很多问题都是好方法，很简单直观（我这里的方法就繁琐很多）。但是忽视了一个重要问题：hull变换后，会丢失信息。显然这就是结果不准确的原因。

我也在answeropencv上进行了咨询，berak给出的comment是 maybe: find principal axes of your shape (PCA) rotate to upright (warp) boundingRect() (image axis aligned) 这个显然也不是很妥当，因为这个结果还需要rotate回去，这会很麻烦。感谢阅读至此，希望有所帮助。

2018年8月29日19:43:01 经过一段时间后反思这个项目，应该说这个算法有一定自创的元素在里面，但是由于这个问题比较小众，所以即使是在answeropencv上面，参与讨论的人也比较少，使得算法有很多不充分的地方在里面：最主要的问题就是在算法的后面部分，多次进行全轮廓循环，使得算法的效率降低。那么，有没有能够提速的方法了？还是之前做的一些数学的研究和在answeropencv上berak的comment提醒了我，注意看：

这里，黑色的是原始的OpenCV的坐标系，红色的是新求出来的坐标系，你花了那么大功夫去算交点，实际上，不如将这个图像旋转为正，将外界矩形算出来，然后再反方向旋转回去。这样的话思路就很清楚了，但是需要花一点修改的时间。但是正是因为这里的思路比较清晰，所以代码写起来，比较流畅，很快我就得到了下面的结果：

1. #include "stdafx.h"

2. #include "opencv2/imgcodecs.hpp"

3. #include "opencv2/highgui.hpp"

4. #include "opencv2/imgproc.hpp"

5. #include "opencv2/photo.hpp"

6. using namespace std;

7. using namespace cv;

8. #define DEBUG FALSE

9. //获得单个点经过旋转后所在精确坐标

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

11.    Point2d preturn;

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

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

14.    return preturn;

15. }

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

17.    Point preturn;

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

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

20.    return preturn;

21. }

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

23. {

24.    //Construct a buffer used by the pca analysis

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

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

27.    {

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

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

30.    }

31.    //Perform PCA analysis

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

33.    //Store the position of the object

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

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

36.    //Store the eigenvalues and eigenvectors

37.    vector<Point2d> eigen_vecs(2);

38.    vector<double> eigen_val(2);

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

40.    {

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

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

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

44.    }

45.    // Draw the principal components

46.    //在轮廓/图像中点绘制小圆

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

48.    计算出直线，在主要方向上绘制直线

49.    //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));

50.    //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));

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

52. }

53. //程序主要部分

54. int main( int argc, char** argv )

55. {

56.    //读入图像，转换为灰度

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

58.    pyrDown(img,img);

59.    pyrDown(img,img);

60.    Mat bw;

61.    bool dRet;

62.    cvtColor(img, bw, COLOR_BGR2GRAY);

63.    //阈值处理

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

65.    //寻找轮廓

66.    vector<vector<Point> > contours;

67.    vector<Vec4i> hierarchy;

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

69.    //轮廓分析，找到

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

71.    {

72.        //计算轮廓大小

73.        double area = contourArea(contours[i]);

74.        //去除过小或者过大的轮廓区域（科学计数法表示）

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

76.        //绘制轮廓

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

78.        //获得轮廓的角度

79.        Point2f* pos = new Point2f();

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

81.        //转换轮廓,并获得极值

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

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

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

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

86.        //将角度转换回去并绘图

87.        Point2f rect_points[4];

88.        rotateRect.points( rect_points );

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

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

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

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

93.        //得出结果    

94.        char cbuf[255];

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

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

97.        sprintf_s(cbuf,"第%d个轮廓,长度%.2f,宽度%.2f像素\n",i,flong,fshort);

98.    }

99.    return 0;

100. }

这段代码中值得一提的是

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

2.    Point preturn;

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

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

5.    return preturn;

6. }

这个函数是直接计算出某一个点在旋转后位置，采用的是数学方法推到，应该算自己创的函数。很多时候，我们并不需要旋转整个图像，而只是要获得图像旋转以后的位置。

反思小结：应该说当时answerOpenCV上就给出了正确的结果提示，但是由于那时我钻在自己的算法里面，没能够接受新的想法；过去一段时间后回顾，才发现了更好的解决方法。但是走弯路并不可怕，只有不断、持续地思考，尽可能将现有的解决方法优化，才可能在面对新的问题的时候有更多的手段、更容易提出创造出“方便书写、效果显著”的算法。此外，基础能力非常重要，如果基础不牢，在创建新算法的时候会遇到更多的困难，毕竟：基础不牢、地动山摇。

link