主成分分析 (PCA) 入門

目次

• 目的
• PCA とは何か?
• 固有ベクトルと固有値はどのように計算されるか?
• ソースコード
• 解説
• 結果

前のチュートリアル: 非線形分離可能データに対するサポートベクターマシン

原著者	Theodore Tsesmelis
互換性	OpenCV >= 3.0

目的

このチュートリアルでは、以下の方法を学ぶ:

• OpenCV のクラス cv::PCA を用いて、物体の向きを算出する。

PCA とは何か?

主成分分析 (PCA) は、データセットの最も重要な特徴を抽出する統計的手法である。

上の図に示すような2次元の点の集合を持っていると考えてほしい。各次元は、関心のある特徴に対応する。ここで、点はランダムな順序で配置されていると主張する人もいるかもしれない。しかしよく見ると、無視しがたい線形パターン(青い線で示される)があることがわかる。PCA の重要な要点は次元削減である。次元削減とは、与えられたデータセットの次元数を減らす処理である。例えば上記の場合、点の集合を1本の線で近似でき、その結果として与えられた点の次元を2次元から1次元へ削減できる。

さらに、点は特徴 1 や特徴 2 の軸方向よりも、青い線に沿った方向に最も大きく分散していることもわかる。これは、青い線に沿った点の位置がわかれば、特徴 1 軸または特徴 2 軸のどこにあるかだけを知っている場合よりも、その点についてより多くの情報を得られることを意味する。

したがって PCA により、データが最も大きく分散する方向を見つけることができる。実際、図中の点の集合に対して PCA を実行した結果は、データセットの主成分である 固有ベクトル (eigenvectors) と呼ばれる2つのベクトルから成る。

各固有ベクトルの大きさは対応する固有値に符号化されており、その主成分に沿ってデータがどれだけ分散しているかを示す。固有ベクトルの始点は、データセット内の全点の中心である。N 次元のデータセットに PCA を適用すると、N 個の N 次元固有ベクトル、N 個の固有値、および1個の N 次元中心点が得られる。理論はこれで十分なので、これらの考え方をどのようにコードへ落とし込むか見ていこう。

固有ベクトルと固有値はどのように計算されるか?

目的は、次元 p の与えられたデータセット X を、より小さい次元 L の別のデータセット Y へ変換することである。言い換えると、行列 X の Karhunen–Loève 変換 (KLT) である行列 Y を求めようとしている。

\[ \mathbf{Y} = \mathbb{K} \mathbb{L} \mathbb{T} \{\mathbf{X}\} \]

データセットを整理する

p 個の変数の観測値の集合から成るデータを持っており、各観測値が L 個の変数のみで記述できるようにデータを削減したいとする(L < p)。さらに、データが n 個のデータベクトル \( x_1...x_n \) の集合として並んでおり、各 \( x_i \) が p 個の変数についての1つのグループ化された観測値を表すとする。

• \( x_1...x_n \) を行ベクトルとして書き、それぞれが p 列を持つようにする。
• それらの行ベクトルを、次元 \( n\times p \) の単一の行列 X にまとめる。

経験平均を計算する

• 各次元 \( j = 1, ..., p \) に沿った経験平均を求める。

• 計算した平均値を、次元 \( p\times 1 \) の経験平均ベクトル u に格納する。

  \[ \mathbf{u[j]} = \frac{1}{n}\sum_{i=1}^{n}\mathbf{X[i,j]} \]

平均からの偏差を計算する

平均の減算は、データを近似する際の平均二乗誤差を最小化する主成分基底を見つける解法に不可欠な部分である。そこで、次のようにデータを中心化していく。

• データ行列 X の各行から、経験平均ベクトル u を減算する。

• 平均を減算したデータを \( n\times p \) 行列 B に格納する。

  \[ \mathbf{B} = \mathbf{X} - \mathbf{h}\mathbf{u^{T}} \]

  ここで h はすべて 1 から成る \( n\times 1 \) の列ベクトルである。

  \[ h[i] = 1, i = 1, ..., n \]

共分散行列を求める

• 行列 B 自身との外積から、\( p\times p \) の経験共分散行列 C を求める。

  \[ \mathbf{C} = \frac{1}{n-1} \mathbf{B^{*}} \cdot \mathbf{B} \]

  ここで * は共役転置演算子である。なお、多くの応用で当てはまるように B が完全に実数から成る場合、「共役転置」は通常の転置と同じになる。

共分散行列の固有ベクトルと固有値を求める

• 共分散行列 C を対角化する固有ベクトルの行列 V を計算する。

  \[ \mathbf{V^{-1}} \mathbf{C} \mathbf{V} = \mathbf{D} \]

  ここで D は C の固有値からなる対角行列である。

• 行列 D は \( p \times p \) の対角行列の形をとる。

  \[ D[k,l] = \left\{\begin{matrix} \lambda_k, k = l \\ 0, k \neq l \end{matrix}\right. \]

  ここで \( \lambda_j \) は共分散行列 C の j 番目の固有値である

• 同じく次元 p x p の行列 V は、それぞれ長さ p の p 個の列ベクトルを含み、これらが共分散行列 C の p 個の固有ベクトルを表す。
• 固有値と固有ベクトルは順序付けられ、対になっている。j 番目の固有値は j 番目の固有ベクトルに対応する。

覚え書き

出典 [1], [2]、および元のチュートリアルを提供してくれた Svetlin Penkov に特に感謝する。

ソースコード

C++

• ダウンロード可能なコード: こちらをクリック
• コードの概要:

   
  #include "opencv2/core.hpp"
  #include "opencv2/imgproc.hpp"
  #include "opencv2/highgui.hpp"
  #include <iostream>
   
  using namespace std;
  using namespace cv;
   
  // 関数宣言
  void drawAxis(Mat&, Point, Point, Scalar, const float);
  double getOrientation(const vector<Point> &, Mat&);
   
  void drawAxis(Mat& img, Point p, Point q, Scalar colour, const float scale = 0.2)
  {
   double angle = atan2( (double) p.y - q.y, (double) p.x - q.x ); // 角度(ラジアン)
   double hypotenuse = sqrt( (double) (p.y - q.y) * (p.y - q.y) + (p.x - q.x) * (p.x - q.x));
   
   // ここで矢印を scale 倍に伸ばす
   q.x = (int) (p.x - scale * hypotenuse * cos(angle));
   q.y = (int) (p.y - scale * hypotenuse * sin(angle));
   line(img, p, q, colour, 1, LINE_AA);
   
   // 矢印の先端(かぎ)を作成
   p.x = (int) (q.x + 9 * cos(angle + CV_PI / 4));
   p.y = (int) (q.y + 9 * sin(angle + CV_PI / 4));
   line(img, p, q, colour, 1, LINE_AA);
   
   p.x = (int) (q.x + 9 * cos(angle - CV_PI / 4));
   p.y = (int) (q.y + 9 * sin(angle - CV_PI / 4));
   line(img, p, q, colour, 1, LINE_AA);
  }
   
  double getOrientation(const vector<Point> &pts, Mat &img)
  {
   // PCA解析で使うバッファを構築
   int sz = static_cast<int>(pts.size());
   Mat data_pts = Mat(sz, 2, CV_64F);
   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;
   }
   
   // PCA解析を実行
   PCA pca_analysis(data_pts, Mat(), PCA::DATA_AS_ROW);
   
   // 物体の中心を格納
   Point cntr = Point(static_cast<int>(pca_analysis.mean.at<double>(0, 0)),
   static_cast<int>(pca_analysis.mean.at<double>(0, 1)));
   
   // 固有値と固有ベクトルを格納
   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);
   }
   
   // 主成分を描画
   circle(img, cntr, 3, Scalar(255, 0, 255), 2);
   Point p1 = cntr + 0.02 * Point(static_cast<int>(eigen_vecs[0].x * eigen_val[0]), static_cast<int>(eigen_vecs[0].y * eigen_val[0]));
   Point p2 = cntr - 0.02 * Point(static_cast<int>(eigen_vecs[1].x * eigen_val[1]), static_cast<int>(eigen_vecs[1].y * eigen_val[1]));
   drawAxis(img, cntr, p1, Scalar(0, 255, 0), 1);
   drawAxis(img, cntr, p2, Scalar(255, 255, 0), 5);
   
   double angle = atan2(eigen_vecs[0].y, eigen_vecs[0].x); // 向き(ラジアン)
   
   return angle;
  }
   
  int main(int argc, char** argv)
  {
   // 画像を読み込む
   CommandLineParser parser(argc, argv, "{@input | pca_test1.jpg | input image}");
   parser.about( "This program demonstrates how to use OpenCV PCA to extract the orientation of an object.\n" );
   parser.printMessage();
   
   Mat src = imread( samples::findFile( parser.get<String>("@input") ) );
   
   // 画像が正しく読み込まれたか確認
   if(src.empty())
   {
   cout << "Problem loading image!!!" << endl;
   return EXIT_FAILURE;
   }
   
   imshow("src", src);
   
   // 画像をグレースケールに変換
   Mat gray;
   cvtColor(src, gray, COLOR_BGR2GRAY);
   
   // 画像を2値画像に変換
   Mat bw;
   threshold(gray, bw, 50, 255, THRESH_BINARY | THRESH_OTSU);
   
   // しきい値処理した画像からすべての輪郭を検出
   vector<vector<Point> > contours;
   findContours(bw, contours, RETR_LIST, CHAIN_APPROX_NONE);
   
   for (size_t i = 0; i < contours.size(); i++)
   {
   // 各輪郭の面積を計算
   double area = contourArea(contours[i]);
   // 小さすぎる、または大きすぎる輪郭は無視
   if (area < 1e2 || 1e5 < area) continue;
   
   // 可視化のためだけに各輪郭を描画
   drawContours(src, contours, static_cast<int>(i), Scalar(0, 0, 255), 2);
   // 各形状の向きを求める
   getOrientation(contours[i], src);
   }
   
   imshow("output", src);
   
   waitKey();
   return EXIT_SUCCESS;
  }

Java

• ダウンロード可能なコード: こちらをクリック
• コードの概要:

  import java.util.ArrayList;
  import java.util.List;
   
  import org.opencv.core.Core;
  import org.opencv.core.CvType;
  import org.opencv.core.Mat;
  import org.opencv.core.MatOfPoint;
  import org.opencv.core.Point;
  import org.opencv.core.Scalar;
  import org.opencv.highgui.HighGui;
  import org.opencv.imgcodecs.Imgcodecs;
  import org.opencv.imgproc.Imgproc;
   
  // このプログラムはOpenCVのPCAを使って物体の向きを抽出する方法を示す。
  class IntroductionToPCA {
   private void drawAxis(Mat img, Point p_, Point q_, Scalar colour, float scale) {
   Point p = new Point(p_.x, p_.y);
   Point q = new Point(q_.x, q_.y);
   double angle = Math.atan2(p.y - q.y, p.x - q.x); // 角度(ラジアン)
   double hypotenuse = Math.sqrt((p.y - q.y) * (p.y - q.y) + (p.x - q.x) * (p.x - q.x));
   
   // ここで矢印を scale 倍に伸ばす
   q.x = (int) (p.x - scale * hypotenuse * Math.cos(angle));
   q.y = (int) (p.y - scale * hypotenuse * Math.sin(angle));
   Imgproc.line(img, p, q, colour, 1, Imgproc.LINE_AA, 0);
   
   // 矢印の先端(かぎ)を作成
   p.x = (int) (q.x + 9 * Math.cos(angle + Math.PI / 4));
   p.y = (int) (q.y + 9 * Math.sin(angle + Math.PI / 4));
   Imgproc.line(img, p, q, colour, 1, Imgproc.LINE_AA, 0);
   
   p.x = (int) (q.x + 9 * Math.cos(angle - Math.PI / 4));
   p.y = (int) (q.y + 9 * Math.sin(angle - Math.PI / 4));
   Imgproc.line(img, p, q, colour, 1, Imgproc.LINE_AA, 0);
   }
   
   private double getOrientation(MatOfPoint ptsMat, Mat img) {
   List<Point> pts = ptsMat.toList();
   // PCA解析で使うバッファを構築
   int sz = pts.size();
   Mat dataPts = new Mat(sz, 2, CvType.CV_64F);
   double[] dataPtsData = new double[(int) (dataPts.total() * dataPts.channels())];
   for (int i = 0; i < dataPts.rows(); i++) {
   dataPtsData[i * dataPts.cols()] = pts.get(i).x;
   dataPtsData[i * dataPts.cols() + 1] = pts.get(i).y;
   }
   dataPts.put(0, 0, dataPtsData);
   
   // PCA解析を実行
   Mat mean = new Mat();
   Mat eigenvectors = new Mat();
   Mat eigenvalues = new Mat();
   Core.PCACompute2(dataPts, mean, eigenvectors, eigenvalues);
   double[] meanData = new double[(int) (mean.total() * mean.channels())];
   mean.get(0, 0, meanData);
   
   // 物体の中心を格納
   Point cntr = new Point(meanData[0], meanData[1]);
   
   // 固有値と固有ベクトルを格納
   double[] eigenvectorsData = new double[(int) (eigenvectors.total() * eigenvectors.channels())];
   double[] eigenvaluesData = new double[(int) (eigenvalues.total() * eigenvalues.channels())];
   eigenvectors.get(0, 0, eigenvectorsData);
   eigenvalues.get(0, 0, eigenvaluesData);
   
   // 主成分を描画
   Imgproc.circle(img, cntr, 3, new Scalar(255, 0, 255), 2);
   Point p1 = new Point(cntr.x + 0.02 * eigenvectorsData[0] * eigenvaluesData[0],
   cntr.y + 0.02 * eigenvectorsData[1] * eigenvaluesData[0]);
   Point p2 = new Point(cntr.x - 0.02 * eigenvectorsData[2] * eigenvaluesData[1],
   cntr.y - 0.02 * eigenvectorsData[3] * eigenvaluesData[1]);
   drawAxis(img, cntr, p1, new Scalar(0, 255, 0), 1);
   drawAxis(img, cntr, p2, new Scalar(255, 255, 0), 5);
   
   double angle = Math.atan2(eigenvectorsData[1], eigenvectorsData[0]); // 向き(ラジアン)
   
   return angle;
   }
   
   public void run(String[] args) {
   // 画像を読み込む
   String filename = args.length > 0 ? args[0] : "../data/pca_test1.jpg";
   Mat src = Imgcodecs.imread(filename);
   
   // 画像が正しく読み込まれたか確認
   if (src.empty()) {
   System.err.println("Cannot read image: " + filename);
   System.exit(0);
   }
   
   Mat srcOriginal = src.clone();
   HighGui.imshow("src", srcOriginal);
   
   // 画像をグレースケールに変換
   Mat gray = new Mat();
   Imgproc.cvtColor(src, gray, Imgproc.COLOR_BGR2GRAY);
   
   // 画像を2値画像に変換
   Mat bw = new Mat();
   Imgproc.threshold(gray, bw, 50, 255, Imgproc.THRESH_BINARY | Imgproc.THRESH_OTSU);
   
   // しきい値処理した画像からすべての輪郭を検出
   List<MatOfPoint> contours = new ArrayList<>();
   Mat hierarchy = new Mat();
   Imgproc.findContours(bw, contours, hierarchy, Imgproc.RETR_LIST, Imgproc.CHAIN_APPROX_NONE);
   
   for (int i = 0; i < contours.size(); i++) {
   // 各輪郭の面積を計算
   double area = Imgproc.contourArea(contours.get(i));
   // 小さすぎる、または大きすぎる輪郭は無視
   if (area < 1e2 || 1e5 < area)
   continue;
   
   // 可視化のためだけに各輪郭を描画
   Imgproc.drawContours(src, contours, i, new Scalar(0, 0, 255), 2);
   // 各形状の向きを求める
   getOrientation(contours.get(i), src);
   }
   
   HighGui.imshow("output", src);
   
   HighGui.waitKey();
   System.exit(0);
   }
  }
   
  public class IntroductionToPCADemo {
   public static void main(String[] args) {
   // OpenCVのネイティブライブラリを読み込む
   System.loadLibrary(Core.NATIVE_LIBRARY_NAME);
   
   new IntroductionToPCA().run(args);
   }
  }

Python

• ダウンロード可能なコード: こちらをクリック
• コードの概要:

  from __future__ import print_function
  from __future__ import division
  import cv2 as cv
  import numpy as np
  import argparse
  from math import atan2, cos, sin, sqrt, pi
   
  def drawAxis(img, p_, q_, colour, scale):
   p = list(p_)
   q = list(q_)
   
   angle = atan2(p[1] - q[1], p[0] - q[0]) # 角度(ラジアン)
   hypotenuse = sqrt((p[1] - q[1]) * (p[1] - q[1]) + (p[0] - q[0]) * (p[0] - q[0]))
   
   # ここで矢印を scale 倍に伸ばす
   q[0] = p[0] - scale * hypotenuse * cos(angle)
   q[1] = p[1] - scale * hypotenuse * sin(angle)
   cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)
   
   # 矢印の先端(かぎ)を作成
   p[0] = q[0] + 9 * cos(angle + pi / 4)
   p[1] = q[1] + 9 * sin(angle + pi / 4)
   cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)
   
   p[0] = q[0] + 9 * cos(angle - pi / 4)
   p[1] = q[1] + 9 * sin(angle - pi / 4)
   cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)
   
   
  def getOrientation(pts, img):
   
   sz = len(pts)
   data_pts = np.empty((sz, 2), dtype=np.float64)
   for i in range(data_pts.shape[0]):
   data_pts[i,0] = pts[i,0,0]
   data_pts[i,1] = pts[i,0,1]
   
   # PCA解析を実行
   mean = np.empty((0))
   mean, eigenvectors, eigenvalues = cv.PCACompute2(data_pts, mean)
   
   # 物体の中心を格納
   cntr = (int(mean[0,0]), int(mean[0,1]))
   
   
   
   cv.circle(img, cntr, 3, (255, 0, 255), 2)
   p1 = (cntr[0] + 0.02 * eigenvectors[0,0] * eigenvalues[0,0], cntr[1] + 0.02 * eigenvectors[0,1] * eigenvalues[0,0])
   p2 = (cntr[0] - 0.02 * eigenvectors[1,0] * eigenvalues[1,0], cntr[1] - 0.02 * eigenvectors[1,1] * eigenvalues[1,0])
   drawAxis(img, cntr, p1, (0, 255, 0), 1)
   drawAxis(img, cntr, p2, (255, 255, 0), 5)
   
   angle = atan2(eigenvectors[0,1], eigenvectors[0,0]) # 向き(ラジアン)
   
   
   return angle
   
   
  parser = argparse.ArgumentParser(description='Code for Introduction to Principal Component Analysis (PCA) tutorial.\
   This program demonstrates how to use OpenCV PCA to extract the orientation of an object.')
  parser.add_argument('--input', help='Path to input image.', default='pca_test1.jpg')
  args = parser.parse_args()
   
  src = cv.imread(cv.samples.findFile(args.input))
  # 画像が正しく読み込まれたか確認
  if src is None:
   print('Could not open or find the image: ', args.input)
   exit(0)
   
  cv.imshow('src', src)
   
  # 画像をグレースケールに変換
  gray = cv.cvtColor(src, cv.COLOR_BGR2GRAY)
   
  # 画像を2値画像に変換
  _, bw = cv.threshold(gray, 50, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)
   
   
   
  contours, _ = cv.findContours(bw, cv.RETR_LIST, cv.CHAIN_APPROX_NONE)
   
  for i, c in enumerate(contours):
   # 各輪郭の面積を計算
   area = cv.contourArea(c)
   # 小さすぎる、または大きすぎる輪郭は無視
   if area < 1e2 or 1e5 < area:
   continue
   
   # 可視化のためだけに各輪郭を描画
   cv.drawContours(src, contours, i, (0, 0, 255), 2)
   # 各形状の向きを求める
   getOrientation(c, src)
   
   
  cv.imshow('output', src)
  cv.waitKey()

覚え書き

PCAを使い、一定の分散を保ちながら次元削減を行う別の例は opencv_source_code/samples/cpp/pca.cpp にある

解説

• 画像を読み込み、2値画像に変換する

ここでは、対象とする物体を検出できるようにするために必要な前処理を行う。

C++

 // Load image
 CommandLineParser parser(argc, argv, "{@input | pca_test1.jpg | input image}");
    parser.about( "This program demonstrates how to use OpenCV PCA to extract the orientation of an object.\n" );
    parser.printMessage();
 
 Mat src = imread( samples::findFile( parser.get<String>("@input") ) );
 
 // Check if image is loaded successfully
 if(src.empty())
    {
        cout << "Problem loading image!!!" << endl;
 return EXIT_FAILURE;
    }
 
 imshow("src", src);
 
 // Convert image to grayscale
 Mat gray;
 cvtColor(src, gray, COLOR_BGR2GRAY);
 
 // Convert image to binary
 Mat bw;
 threshold(gray, bw, 50, 255, THRESH_BINARY | THRESH_OTSU);

Java

 // Load image
 String filename = args.length > 0 ? args[0] : "../data/pca_test1.jpg";
        Mat src = Imgcodecs.imread(filename);
 
 // Check if image is loaded successfully
 if (src.empty()) {
            System.err.println("Cannot read image: " + filename);
            System.exit(0);
        }
 
        Mat srcOriginal = src.clone();
        HighGui.imshow("src", srcOriginal);
 
 // Convert image to grayscale
        Mat gray = new Mat();
        Imgproc.cvtColor(src, gray, Imgproc.COLOR_BGR2GRAY);
 
 // Convert image to binary
        Mat bw = new Mat();
        Imgproc.threshold(gray, bw, 50, 255, Imgproc.THRESH_BINARY | Imgproc.THRESH_OTSU);

Python

parser = argparse.ArgumentParser(description='Code for Introduction to Principal Component Analysis (PCA) tutorial.\
                                              This program demonstrates how to use OpenCV PCA to extract the orientation of an object.')
parser.add_argument('--input', help='Path to input image.', default='pca_test1.jpg')
args = parser.parse_args()
 
src = cv.imread(cv.samples.findFile(args.input))
# Check if image is loaded successfully
if src is None:
    print('Could not open or find the image: ', args.input)
    exit(0)
 
cv.imshow('src', src)
 
# Convert image to grayscale
gray = cv.cvtColor(src, cv.COLOR_BGR2GRAY)
 
# Convert image to binary
_, bw = cv.threshold(gray, 50, 255, cv.THRESH_BINARY | cv.THRESH_OTSU)

• 対象とする物体を抽出する

次に、サイズによって輪郭を検出・フィルタリングし、残った輪郭の向きを求める。

C++

 // Find all the contours in the thresholded image
    vector<vector<Point> > contours;
 findContours(bw, contours, RETR_LIST, CHAIN_APPROX_NONE);
 
 for (size_t i = 0; i < contours.size(); i++)
    {
 // Calculate the area of each contour
 double area = contourArea(contours[i]);
 // Ignore contours that are too small or too large
 if (area < 1e2 || 1e5 < area) continue;
 
 // Draw each contour only for visualisation purposes
 drawContours(src, contours, static_cast<int>(i), Scalar(0, 0, 255), 2);
 // Find the orientation of each shape
        getOrientation(contours[i], src);
    }

Java

 // Find all the contours in the thresholded image
        List<MatOfPoint> contours = new ArrayList<>();
        Mat hierarchy = new Mat();
        Imgproc.findContours(bw, contours, hierarchy, Imgproc.RETR_LIST, Imgproc.CHAIN_APPROX_NONE);
 
 for (int i = 0; i < contours.size(); i++) {
 // Calculate the area of each contour
 double area = Imgproc.contourArea(contours.get(i));
 // Ignore contours that are too small or too large
 if (area < 1e2 || 1e5 < area)
 continue;
 
 // Draw each contour only for visualisation purposes
            Imgproc.drawContours(src, contours, i, new Scalar(0, 0, 255), 2);
 // Find the orientation of each shape
            getOrientation(contours.get(i), src);
        }

Python

contours, _ = cv.findContours(bw, cv.RETR_LIST, cv.CHAIN_APPROX_NONE)
 
for i, c in enumerate(contours):
 # Calculate the area of each contour
    area = cv.contourArea(c)
 # Ignore contours that are too small or too large
 if area < 1e2 or 1e5 < area:
 continue
 
 # Draw each contour only for visualisation purposes
 cv.drawContours(src, contours, i, (0, 0, 255), 2)
 # Find the orientation of each shape
    getOrientation(c, src)

• 向きを抽出する

向きは getOrientation() 関数の呼び出しによって抽出され、この関数がPCAの処理全体を行う。

C++

 //Construct a buffer used by the pca analysis
 int sz = static_cast<int>(pts.size());
 Mat data_pts = Mat(sz, 2, CV_64F);
 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(), PCA::DATA_AS_ROW);
 
 //Store the center of the object
 Point cntr = Point(static_cast<int>(pca_analysis.mean.at<double>(0, 0)),
 static_cast<int>(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);
    }

Java

 // Construct a buffer used by the pca analysis
 int sz = pts.size();
        Mat dataPts = new Mat(sz, 2, CvType.CV_64F);
 double[] dataPtsData = new double[(int) (dataPts.total() * dataPts.channels())];
 for (int i = 0; i < dataPts.rows(); i++) {
            dataPtsData[i * dataPts.cols()] = pts.get(i).x;
            dataPtsData[i * dataPts.cols() + 1] = pts.get(i).y;
        }
        dataPts.put(0, 0, dataPtsData);
 
 // Perform PCA analysis
        Mat mean = new Mat();
        Mat eigenvectors = new Mat();
        Mat eigenvalues = new Mat();
        Core.PCACompute2(dataPts, mean, eigenvectors, eigenvalues);
 double[] meanData = new double[(int) (mean.total() * mean.channels())];
 mean.get(0, 0, meanData);
 
 // Store the center of the object
 Point cntr = new Point(meanData[0], meanData[1]);
 
 // Store the eigenvalues and eigenvectors
 double[] eigenvectorsData = new double[(int) (eigenvectors.total() * eigenvectors.channels())];
 double[] eigenvaluesData = new double[(int) (eigenvalues.total() * eigenvalues.channels())];
        eigenvectors.get(0, 0, eigenvectorsData);
        eigenvalues.get(0, 0, eigenvaluesData);

Python

    sz = len(pts)
    data_pts = np.empty((sz, 2), dtype=np.float64)
 for i in range(data_pts.shape[0]):
        data_pts[i,0] = pts[i,0,0]
        data_pts[i,1] = pts[i,0,1]
 
 # Perform PCA analysis
    mean = np.empty((0))
    mean, eigenvectors, eigenvalues = cv.PCACompute2(data_pts, mean)
 
 # Store the center of the object
    cntr = (int(mean[0,0]), int(mean[0,1]))

まず、データをサイズ n x 2 の行列に並べる必要がある。ここで n はデータ点の数である。その後、PCA解析を行うことができる。計算された平均(すなわち重心)は cntr 変数に格納され、固有ベクトルと固有値はそれぞれ対応する std::vector に格納される。

• 結果を可視化する

最終的な結果は drawAxis() 関数によって可視化される。この関数では主成分が直線で描画され、各固有ベクトルはその固有値で乗算され、平均の位置へ平行移動される。

C++

 // Draw the principal components
 circle(img, cntr, 3, Scalar(255, 0, 255), 2);
 Point p1 = cntr + 0.02 * Point(static_cast<int>(eigen_vecs[0].x * eigen_val[0]), static_cast<int>(eigen_vecs[0].y * eigen_val[0]));
 Point p2 = cntr - 0.02 * Point(static_cast<int>(eigen_vecs[1].x * eigen_val[1]), static_cast<int>(eigen_vecs[1].y * eigen_val[1]));
    drawAxis(img, cntr, p1, Scalar(0, 255, 0), 1);
    drawAxis(img, cntr, p2, Scalar(255, 255, 0), 5);
 
 double angle = atan2(eigen_vecs[0].y, eigen_vecs[0].x); // orientation in radians

Java

 // Draw the principal components
        Imgproc.circle(img, cntr, 3, new Scalar(255, 0, 255), 2);
 Point p1 = new Point(cntr.x + 0.02 * eigenvectorsData[0] * eigenvaluesData[0],
                cntr.y + 0.02 * eigenvectorsData[1] * eigenvaluesData[0]);
 Point p2 = new Point(cntr.x - 0.02 * eigenvectorsData[2] * eigenvaluesData[1],
                cntr.y - 0.02 * eigenvectorsData[3] * eigenvaluesData[1]);
        drawAxis(img, cntr, p1, new Scalar(0, 255, 0), 1);
        drawAxis(img, cntr, p2, new Scalar(255, 255, 0), 5);
 
 double angle = Math.atan2(eigenvectorsData[1], eigenvectorsData[0]); // orientation in radians

Python

 cv.circle(img, cntr, 3, (255, 0, 255), 2)
    p1 = (cntr[0] + 0.02 * eigenvectors[0,0] * eigenvalues[0,0], cntr[1] + 0.02 * eigenvectors[0,1] * eigenvalues[0,0])
    p2 = (cntr[0] - 0.02 * eigenvectors[1,0] * eigenvalues[1,0], cntr[1] - 0.02 * eigenvectors[1,1] * eigenvalues[1,0])
    drawAxis(img, cntr, p1, (0, 255, 0), 1)
    drawAxis(img, cntr, p2, (255, 255, 0), 5)
 
    angle = atan2(eigenvectors[0,1], eigenvectors[0,0]) # orientation in radians

C++

 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));
 line(img, p, q, colour, 1, 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));
 line(img, p, q, colour, 1, LINE_AA);
 
    p.x = (int) (q.x + 9 * cos(angle - CV_PI / 4));
    p.y = (int) (q.y + 9 * sin(angle - CV_PI / 4));
 line(img, p, q, colour, 1, LINE_AA);

Java

 double angle = Math.atan2(p.y - q.y, p.x - q.x); // angle in radians
 double hypotenuse = Math.sqrt((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 * Math.cos(angle));
        q.y = (int) (p.y - scale * hypotenuse * Math.sin(angle));
        Imgproc.line(img, p, q, colour, 1, Imgproc.LINE_AA, 0);
 
 // create the arrow hooks
        p.x = (int) (q.x + 9 * Math.cos(angle + Math.PI / 4));
        p.y = (int) (q.y + 9 * Math.sin(angle + Math.PI / 4));
        Imgproc.line(img, p, q, colour, 1, Imgproc.LINE_AA, 0);
 
        p.x = (int) (q.x + 9 * Math.cos(angle - Math.PI / 4));
        p.y = (int) (q.y + 9 * Math.sin(angle - Math.PI / 4));
        Imgproc.line(img, p, q, colour, 1, Imgproc.LINE_AA, 0);

Python

    angle = atan2(p[1] - q[1], p[0] - q[0]) # angle in radians
    hypotenuse = sqrt((p[1] - q[1]) * (p[1] - q[1]) + (p[0] - q[0]) * (p[0] - q[0]))
 
 # Here we lengthen the arrow by a factor of scale
    q[0] = p[0] - scale * hypotenuse * cos(angle)
    q[1] = p[1] - scale * hypotenuse * sin(angle)
 cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)
 
 # create the arrow hooks
    p[0] = q[0] + 9 * cos(angle + pi / 4)
    p[1] = q[1] + 9 * sin(angle + pi / 4)
 cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)
 
    p[0] = q[0] + 9 * cos(angle - pi / 4)
    p[1] = q[1] + 9 * sin(angle - pi / 4)
 cv.line(img, (int(p[0]), int(p[1])), (int(q[0]), int(q[1])), colour, 1, cv.LINE_AA)

結果

このコードは画像を開き、検出した対象物体の向きを求め、検出した対象物体の輪郭、中心点、および抽出した向きに対応するx軸とy軸を描画して結果を可視化する。