主成分分析 PCA (Principal Component Analysis) を利用すると、多次元のデータから興味のある特徴のみを抽出して、より低次元のデータに変換できる可能性があります。PCA の応用例の一つとして、物体の方向検出があります。画像をグレースケールに変換して物体の輪郭を検出した後に、輪郭を形成する各点を x座標、y座標を持つ二次元データとみなします。この二次元データに対して PCA を適用して、一次元の情報に圧縮することで方向を検出します。

物体の方向検出

OpenCV に付属のグレースケール画像を利用して、物体の方向を検出してみます。

#include <opencv2/opencv.hpp>
#include <iostream>

void drawAxis(cv::Mat& img, cv::Point p, cv::Point q, cv::Scalar colour, const float scale = 0.2) {
    double angle = atan2( (double) p.y - q.y, (double) p.x - q.x ); // angle in radians
    double hypotenuse = sqrt( (double) (p.y - q.y) * (p.y - q.y) + (p.x - q.x) * (p.x - q.x));
    // Here we lengthen the arrow by a factor of scale
    q.x = (int) (p.x - scale * hypotenuse * cos(angle));
    q.y = (int) (p.y - scale * hypotenuse * sin(angle));
    cv::line(img, p, q, colour, 1, cv::LINE_AA);
    // create the arrow hooks
    p.x = (int) (q.x + 9 * cos(angle + CV_PI / 4));
    p.y = (int) (q.y + 9 * sin(angle + CV_PI / 4));
    cv::line(img, p, q, colour, 1, cv::LINE_AA);
    p.x = (int) (q.x + 9 * cos(angle - CV_PI / 4));
    p.y = (int) (q.y + 9 * sin(angle - CV_PI / 4));
    cv::line(img, p, q, colour, 1, cv::LINE_AA);
}

double getOrientation(const std::vector<cv::Point> &pts, cv::Mat &img) {
    // 主成分分析のために、対象となる輪郭のみを別のバッファにコピーします。
    cv::Mat data_pts = cv::Mat(pts.size(), 2, CV_64F); // [pts.size() x 2] 行列
    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;
    }

    // 主成分分析の実行
    cv::PCA pca_analysis(data_pts, cv::Mat(), cv::PCA::DATA_AS_ROW);

    // [分析結果] 中心の座標
    cv::Point cntr = cv::Point(pca_analysis.mean.at<double>(0, 0),
                               pca_analysis.mean.at<double>(0, 1));
    std::cout << pca_analysis.mean << std::endl;

    // [分析結果] 固有値と固有ベクトル
    std::vector<cv::Point2d> eigen_vecs(2);
    std::vector<double> eigen_val(2);
    for (int i = 0; i < 2; i++) {
        eigen_vecs[i] = cv::Point2d(pca_analysis.eigenvectors.at<double>(i, 0),
                                    pca_analysis.eigenvectors.at<double>(i, 1));
        eigen_val[i] = pca_analysis.eigenvalues.at<double>(i);
    }
    std::cout << pca_analysis.eigenvectors << std::endl;
    std::cout << pca_analysis.eigenvalues << std::endl; // 固有値が小さい側について、次元の削除を検討します。

    // 分析結果を可視化のために描画します。
    cv::circle(img, cntr, 3, cv::Scalar(255, 0, 255), 2);
    cv::Point p1 = cntr + 0.02 * cv::Point(eigen_vecs[0].x * eigen_val[0], eigen_vecs[0].y * eigen_val[0]);
    cv::Point p2 = cntr - 0.02 * cv::Point(eigen_vecs[1].x * eigen_val[1], eigen_vecs[1].y * eigen_val[1]);
    drawAxis(img, cntr, p1, cv::Scalar(0, 255, 0), 1);
    drawAxis(img, cntr, p2, cv::Scalar(255, 255, 0), 5);

    // 物体の方向をラジアンで返します。
    double angle = atan2(eigen_vecs[0].y, eigen_vecs[0].x);
    return angle;
}

int main() {
    cv::Mat src = cv::imread("pca_test1.jpg", -1);
    if(src.empty()) {
        return -1;
    }

    // グレースケール画像に変換
    cv::Mat gray;
    cv::cvtColor(src, gray, cv::COLOR_BGR2GRAY);

    // 50 を閾値として黒 0 か白 255 に二値化します。
    cv::Mat bw;
    cv::threshold(gray, bw, 50, 255, cv::THRESH_BINARY | cv::THRESH_OTSU);
    cv::imshow("gray", gray);
    cv::imshow("bw", bw);

    // 輪郭を検出します。
    std::vector<std::vector<cv::Point> > contours;
    cv::findContours(bw, contours, cv::RETR_LIST, cv::CHAIN_APPROX_NONE);

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

        // 面積が小さすぎる場合と多きすぎる場合を除外します。
        double area = cv::contourArea(contours[i]);
        if (area < 1e2 || 1e5 < area) {
            continue;
        }
        // 輪郭を描画します。
        cv::drawContours(src, contours, i, cv::Scalar(0, 0, 255), 2);
        // 各輪郭について、それぞれ方向を検出します。
        getOrientation(contours[i], src);
    }

    cv::imshow("output", src);
    cv::waitKey(0);
    return 0;
}

各輪郭について、中心点の座標、固有ベクトル、固有値

[430.2620865139949, 407.7417302798983]
[0.966420963693683, -0.256964045993545;
 0.256964045993545, 0.966420963693683]
[13507.35794872479;
 499.8921548692795]

[439.8099489795918, 326.015306122449]
[0.9724910769422836, -0.232940132368034;
 0.232940132368034, 0.9724910769422836]
[13275.51140165105;
 486.8132140689587]

[433.4980595084089, 239.7697283311772]
[0.9753789596742241, -0.2205354507212575;
 0.2205354507212575, 0.9753789596742241]
[12837.25434334891;
 463.1703121908179]

[420.184655396619, 169.1131339401821]
[0.985849803253799, -0.1676310395614297;
 0.1676310395614297, 0.985849803253799]
[12496.98682699333;
 430.1275244772632]

[191.3530183727034, 291.990813648294]
[0.06092594794314996, 0.998142288888327;
 0.998142288888327, -0.06092594794314996]
[11958.1181863091;
 463.3791545744922]

[407.8350923482849, 90.51187335092348]
[0.9854469086343497, -0.1699834999728027;
 0.1699834999728027, 0.9854469086343497]
[12096.51049439792;
 414.2359167938764]