OpenCV
开源计算机视觉库
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
主成分分析 (PCA) 入门

上一教程: 用于非线性可分离数据的支持向量机

原作者Theodore Tsesmelis
兼容性OpenCV >= 3.0

目标

在本教程中,您将学习如何

  • 使用 OpenCV 类 cv::PCA 计算物体的方向。

什么是 PCA?

主成分分析 (PCA) 是一种统计程序,用于提取数据集最重要的特征。

假设您有一组 2D 点,如上图所示。每个维度对应于您感兴趣的一个特征。有人可能会认为这些点是随机排列的。但是,如果您仔细观察,您会发现存在一个线性模式(由蓝线指示),很难忽视。PCA 的一个关键点是降维。降维是减少给定数据集维度数量的过程。例如,在上述情况下,可以将点集近似为一条直线,从而将给定点的维度从 2D 降到 1D。

此外,您还可以看到,沿着蓝线方向,点的变化比沿着特征 1 或特征 2 轴方向的变化更大。这意味着,如果您知道点沿蓝线的位移,那么您比只知道它在特征 1 轴或特征 2 轴上的位置时,拥有更多关于该点的信息。

因此,PCA 允许我们找到数据变化最大的方向。事实上,对图中点集运行 PCA 的结果包括 2 个称为 *特征向量* 的向量,它们是数据集的 *主成分*。

每个特征向量的长度由相应的特征值编码,并指示数据沿主成分变化的程度。特征向量的起点是数据集中所有点的中心。将 PCA 应用于 N 维数据集将产生 N 个 N 维特征向量、N 个特征值和 1 个 N 维中心点。理论讲得够多了,让我们看看如何将这些想法转化为代码。

如何计算特征向量和特征值?

目标是将给定维度为 *p* 的数据集 **X** 变换为另一个维度较小的数据集 **Y**,维度为 *L*。等效地,我们正在寻找矩阵 **Y**,其中 **Y** 是矩阵 **X** 的 *Karhunen–Loève 变换* (KLT)

Y=KLT{X}

组织数据集

假设您拥有的数据包含对 *p* 个变量的一组观测值,并且您希望减少数据,以便每个观测值都可以仅用 *L* 个变量来描述,其中 *L* < *p*。进一步假设,数据排列为一组 *n* 个数据向量 x1...xn,每个 xi 代表对 *p* 个变量的单个分组观测值。

  • x1...xn 写成行向量,每个行向量都有 *p* 列。
  • 将行向量放入一个维度为 n×p 的单个矩阵 **X** 中。

计算经验均值

  • 沿每个维度 j=1,...,p 找到经验均值。
  • 将计算出的均值放入维度为 p×1 的经验均值向量 **u** 中。

    u[j]=1ni=1nX[i,j]

计算与均值的偏差

均值减法是寻找最小化数据近似均方误差的主成分基的解决方案中不可或缺的一部分。因此,我们按如下方式对数据进行中心化:

  • 从数据矩阵 **X** 的每一行减去经验均值向量 **u**。
  • 将均值减去后的数据存储在 n×p 矩阵 **B** 中。

    B=XhuT

    其中 **h** 是一个 n×1 的全为 1 的列向量

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

找到协方差矩阵

  • 从矩阵 **B** 与自身的外部积中找到 p×p 的经验协方差矩阵 **C**

    C=1n1BB

    其中 * 是共轭转置运算符。请注意,如果 B 完全由实数组成,这在许多应用中都是如此,“共轭转置”与常规转置相同。

找到协方差矩阵的特征向量和特征值

  • 计算特征向量矩阵 **V**,它对协方差矩阵 **C** 进行对角化

    V1CV=D

    其中 **D** 是 **C** 的特征值的对角矩阵。

  • 矩阵 **D** 将采用 p×p 对角矩阵的形式

    D[k,l]={λk,k=l0,kl

    这里,λj 是协方差矩阵 **C** 的第 *j* 个特征值

  • 矩阵 **V** 的维度也是 *p* x *p*,包含 *p* 个列向量,每个向量的长度为 *p*,它们代表协方差矩阵 **C** 的 *p* 个特征向量。
  • 特征值和特征向量是有序且成对的。第 *j* 个特征值对应于第 *j* 个特征向量。
注意
来源 [1][2],并特别感谢 Svetlin Penkov 提供的原始教程。

源代码

  • 可下载代码: 点击 此处
  • 代码概览
    #include "opencv2/core.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( "本程序演示如何使用OpenCV PCA提取物体的方向。\n" );
    parser.printMessage();
    Mat src = imread( samples::findFile( parser.get<String>("@input") ) );
    // 检查图像是否加载成功
    if(src.empty())
    {
    cout << "图像加载失败!!!" << endl;
    return EXIT_FAILURE;
    }
    imshow("src", src);
    // 将图像转换为灰度图像
    Mat gray;
    cvtColor(src, gray, COLOR_BGR2GRAY);
    // 将图像转换为二值图像
    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);
    return EXIT_SUCCESS;
    }
    用于命令行解析。
    定义 utility.hpp:890
    n维密集数组类
    定义 mat.hpp:829
    _Tp & at(int i0=0)
    返回对指定数组元素的引用。
    bool empty() const
    如果数组没有元素,则返回 true。
    int rows
    行和列的数量,当矩阵维度超过2维时为(-1, -1)
    定义 mat.hpp:2155
    主成分分析。
    定义 core.hpp:2498
    _Tp y
    点的y坐标
    定义 types.hpp:202
    _Tp x
    点的x坐标
    定义 types.hpp:201
    std::string String
    定义 cvstd.hpp:151
    #define CV_64F
    定义 interface.h:79
    #define CV_PI
    定义 cvdef.h:380
    void imshow(const String &winname, InputArray mat)
    在指定的窗口中显示图像。
    int waitKey(int delay=0)
    等待按键按下。
    CV_EXPORTS_W Mat imread(const String &filename, int flags=IMREAD_COLOR_BGR)
    从文件中加载图像。
    void cvtColor(InputArray src, OutputArray dst, int code, int dstCn=0, AlgorithmHint hint=cv::ALGO_HINT_DEFAULT)
    将图像从一个颜色空间转换为另一个颜色空间。
    void drawContours(InputOutputArray image, InputArrayOfArrays contours, int contourIdx, const Scalar &color, int thickness=1, int lineType=LINE_8, InputArray hierarchy=noArray(), int maxLevel=INT_MAX, Point offset=Point())
    绘制轮廓轮廓或填充轮廓。
    double threshold(InputArray src, OutputArray dst, double thresh, double maxval, int type)
    对每个数组元素应用固定级别的阈值。
    double contourArea(InputArray contour, bool oriented=false)
    计算轮廓面积。
    void findContours(InputArray image, OutputArrayOfArrays contours, OutputArray hierarchy, int mode, int method, Point offset=Point())
    在二值图像中查找轮廓。
    int main(int argc, char *argv[])
    定义 highgui_qt.cpp:3
    定义 core.hpp:107
    STL 命名空间。
注意
另一个使用PCA进行降维(同时保持一定方差)的例子可以在 opencv_source_code/samples/cpp/pca.cpp 找到。

解释

  • 读取图像并将其转换为二值图像

这里我们应用必要的预处理步骤,以便能够检测感兴趣的目标。

// 加载图像
CommandLineParser parser(argc, argv, "{@input | pca_test1.jpg | input image}");
parser.about( "本程序演示如何使用OpenCV PCA提取物体的方向。\n" );
parser.printMessage();
Mat src = imread( samples::findFile( parser.get<String>("@input") ) );
// 检查图像是否加载成功
if(src.empty())
{
cout << "图像加载失败!!!" << endl;
return EXIT_FAILURE;
}
imshow("src", src);
// 将图像转换为灰度图像
Mat gray;
cvtColor(src, gray, COLOR_BGR2GRAY);
// 将图像转换为二值图像
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);
}
  • 提取方向

方向通过调用getOrientation()函数提取,该函数执行所有PCA过程。

// 构造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);
}

首先,需要将数据排列成一个n x 2的矩阵,其中n是我们拥有的数据点的数量。然后我们可以执行PCA分析。计算出的均值(即质心)存储在cntr变量中,特征向量和特征值存储在相应的std::vector中。

  • 可视化结果

最终结果通过drawAxis()函数可视化,其中主成分以线的方式绘制,每个特征向量乘以其特征值并平移到均值位置。

// 绘制主成分
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); // 方向(弧度)
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);

结果

代码打开图像,查找检测到的目标对象的方向,然后通过绘制检测到的目标对象的轮廓、中心点以及关于提取方向的x轴和y轴来可视化结果。