#include<opencv2\opencv.hpp>
using namespace cv;
using namespace std;

double calcPCAorientation(vector<Point>&pts, Mat &image);

int main(int arc, char** argv)
{
    Mat src = imread("3.jpg");
    namedWindow("input", WINDOW_AUTOSIZE);
    imshow("input", src);
    Mat gray, binary;
    cvtColor(src, gray, COLOR_BGR2GRAY);
    threshold(gray, binary, 0, 255, THRESH_BINARY | THRESH_OTSU);
    //imshow("binary", binary);

    vector<vector<Point>>contours;
    findContours(binary, contours, RETR_LIST, CHAIN_APPROX_NONE);
    Mat result = src.clone();
    for (int i = 0; i < contours.size(); i++) {
        double area = contourArea(contours[i]);
        if (area > 1e5 || area < 1e2)continue;
        drawContours(result, contours, i, Scalar(0, 0, 255),1);
        calcPCAorientation(contours[i], result);
    }
    imshow("result", result);

    waitKey(0);
    return 0;
}

double calcPCAorientation(vector<Point>&pts, Mat &image) 
{
    int size = pts.size();
    Mat data_pts = Mat(size, 2, CV_64FC1);
    for (int i = 0; i < size; i++) {
        data_pts.at<double>(i, 0) = pts[i].x;
        data_pts.at<double>(i, 1) = pts[i].y;
    }
    //PCA
    PCA pca_analysis(data_pts, Mat(), PCA::DATA_AS_ROW);
    Point cnt = Point(pca_analysis.mean.at<double>(0, 0), pca_analysis.mean.at<double>(0, 1));//平均值
    circle(image, cnt, 2, Scalar(0, 255, 0), 2);

    vector<Point2d>vecs(2);//特征向量
    vector<double>vals(2);//特征值
    for (int i = 0; i < 2; i++) 
    {
        vals[i] = pca_analysis.eigenvalues.at<double>(i, 0);
        vecs[i] = Point2d(pca_analysis.eigenvectors.at<double>(i, 0), pca_analysis.eigenvectors.at<double>(i, 1));
    }
    Point p1 = cnt + Point(static_cast<int>(vecs[0].x * 320), static_cast<int>(vecs[0].y * 315));//原本是300
    Point p2 = cnt - Point(static_cast<int>(vecs[1].x * 100), static_cast<int>(vecs[1].y * 100));//50
    //即将原直角坐标系中的数据在下面新的一组基中表示，使得数据特征最大化
    line(image, cnt, p1, Scalar(255, 0, 0), 1);//主特征方向
    line(image, cnt, p2, Scalar(0, 255, 0), 1);//与主特征垂直方向
    double angle = atan2(vecs[0].y, vecs[0].x);
    printf("angle:%.2f\n", 180 * angle / CV_PI);
    printf("%d,%d\n", vals[0], vals[1]);
    return angle;
}