OpenCV  
开源计算机视觉库
加载中...
搜索中...
无匹配结果
主成分分析 (PCA) 简介

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

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

目标

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

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

什么是 PCA?

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

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

此外,您还可以看到,与沿着特征 1 或特征 2 轴的变化相比,点沿着蓝线变化更大。这意味着如果您知道点沿着蓝线的坐标,那么您将获得比仅知道它在特征 1 轴或特征 2 轴上的位置更多的信息。

因此,PCA 使我们能够找到数据变化最大的方向。事实上,在点集图中运行 PCA 的结果包含两个称为特征向量的向量,它们是数据集的主成分

每个特征向量的尺寸在相应的特征值中编码,并表示数据沿着主成分变化的程度。特征向量的起点是数据集中所有点的中心。将 PCA 应用于 N 维数据集会产生 N 个 N 维特征向量、N 个特征值和 1 个 N 维中心点。理论足够了,让我们看看如何将这些想法应用到代码中。

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

目标是将给定维度为p的数据集X转换为维度为L的替代数据集Y。等效地,我们正在寻找矩阵Y,其中Y是矩阵X卡尔洪-洛夫变换 (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是一个所有元素均为 1 的n×1列向量

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

找到协方差矩阵

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

    C=1n1BB

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

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

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

    V1CV=D

    其中DC的特征值的对角矩阵。

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

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

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

  • 矩阵V,也为p x p维,包含p个列向量,每个向量长度为p,它们代表协方差矩阵Cp个特征向量。
  • 特征值和特征向量是有序且配对的。第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 | 输入图像}");
    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:820
    n 维密集数组类
    定义 mat.hpp:812
    _Tp & at(int i0=0)
    返回对指定数组元素的引用。
    bool empty() const
    如果数组没有元素,则返回 true。
    int rows
    行数和列数,或者当矩阵具有超过 2 个维度时为 (-1, -1)
    定义 mat.hpp:2138
    主成分分析。
    定义 core.hpp:2513
    _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)
    从文件加载图像。
    void cvtColor(InputArray src, OutputArray dst, int code, int dstCn=0)
    将图像从一个颜色空间转换为另一个颜色空间。
    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:102
    STL 命名空间。
注意
另一个使用 PCA 进行降维,同时保持一定方差的例子可以在 opencv_source_code/samples/cpp/pca.cpp 中找到。

解释

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

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

// 加载图像
CommandLineParser parser(argc, argv, "{@input | pca_test1.jpg | 输入图像}");
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 轴来可视化结果。