PCA的算法应该是比较简单的，一个基于方差的投影降维方式，降维的维度采用方差最大的方式（方差大小反应样本的离散程度，越离散可分性理论上应该越好），PCA降维可以提高分类器的分类效果。本主题使用OpenCV提供的例子，说明PCA可以判别物体的方向。本主题包含两个例子的实现代码：
  1. 鸢尾花的2个特征数据集的分散方向可视化；
  2. 物体方向判别可视化；

---

提示：OpenCV中PCA不在ml模块，而是在core模块中。

---

一.鸢尾花的数据集分散方向可视化

效果

鸢尾花第1与3个特征的离散方向可视化效果

代码

/*
    这个例子中，高宽的比例最好为1，因为高宽比例不一致，主成分的向量的正交，在变换后，会变成不正交的效果
 */

#include <opencv2/opencv.hpp>
#include <iostream>
#include <cmath>
#include <fstream>
#include <sstream>
// 函数定义
void drawAxis(cv::Mat&, cv::Point, cv::Point, cv::Scalar, const float);
// 加载鸢尾花数据集
void loadIrisFromTxt(cv::String, cv::Mat&, cv::Mat&);
void visual_pca(cv::Mat,cv::Mat&, int=512,int h=512);

// 主程序
int main(int argc, const char** argv) {
    // 读取鸢尾花数据集
    cv::Mat data,label;
    loadIrisFromTxt("../iris.txt", data, label);
    // std::cout << data << std::endl;
    // 取其中两列处理
    cv::Mat newData;
    cv::hconcat(data.rowRange(50,100).col(0), data.rowRange(50,100).col(2), newData);

    // 用来缓存所有的图像
    cv::Mat img_out;
    visual_pca(newData,img_out, 600, 600);   // 如果高宽不一致，会产生仿射变换
    cv::imwrite("../iris.png",img_out);
    return 0;
}

void drawAxis(cv::Mat& img, cv::Point p, cv::Point q, cv::Scalar colour, const float scale = 0.2){
    /*
        cv::Mat& img        ：需要绘制的图像
        cv::Point p         ：绘制坐标轴线段的两个点之一：坐标轴的起点
        cv::Point q         ：绘制坐标轴线段的两个点之二：坐标轴的终点
        cv::Scalar colour   ：绘制的颜色
        const float scale   ：坐标轴的放大因子
     */

    // 计算两个点与原点形成的向量的夹角（用于延长线段的计算）
    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延长,并绘制（p不动，按照原来pq的联线方向，增加线段长度）
    q.x = (int) (p.x - scale * hypotenuse * cos(angle));
    q.y = (int) (p.y - scale * hypotenuse * sin(angle));
    line(img, p, q, colour, 1, cv::LINE_AA);

    // 绘制坐标轴的箭头（以终点q为起点的两个小线段，与原来线段成45°夹角），箭头线段9像素
    int len = 9;
    // 箭头：第一个线段(逆时针旋转45°)
    p.x = (int) (q.x + len * cos(angle + CV_PI / 4));
    p.y = (int) (q.y + len * sin(angle + CV_PI / 4));
    line(img, p, q, colour, 1, cv::LINE_AA);
    // 箭头：第二个线段(顺时针旋转45°)
    p.x = (int) (q.x + len * cos(angle - CV_PI / 4));
    p.y = (int) (q.y + len * sin(angle - CV_PI / 4));
    line(img, p, q, colour, 1, cv::LINE_AA);
}


// 从txt文件加载鸢尾花数据
void loadIrisFromTxt(cv::String filename, cv::Mat &data, cv::Mat &label){
    int rows = 150;
    int cols = 4;
    data.create(rows, cols, CV_32FC1);
    label.create(rows, 1, CV_32SC1);
    // 打开文件
    std::ifstream ifs(filename, std::ifstream::in);
    char buf[256];
    int row_id = 0;
    while(ifs.good()){
        bzero(buf, sizeof(buf));
        ifs.getline(buf, sizeof(buf)-1);
        std::istringstream istr(buf);
        float a, b, c, d;
        char buffer[256] = {0};
        char _;
        istr >> a >> _ >> b >> _ >> c >> _ >> d >> _ >> buffer;
        // std::cout << a << "," << b << "," << c << "," << d << std::endl;
        cv::Mat row_one = (cv::Mat_<float>(1,4) << a, b, c, d);
        row_one.copyTo(data.row(row_id));
        
        if(strcmp(buffer, "Iris-setosa") == 0){
            label.at<int>(row_id, 0) = 0;
        }
        else if (strcmp(buffer, "Iris-versicolor") == 0){
            label.at<int>(row_id, 0) = 1;
        }
        else if (! strcmp(buffer, "Iris-virginica")){
            label.at<int>(row_id, 0) = 2;
        }
        row_id ++;
    }
    // 关闭文件
    ifs.close();
};

// 可视化分类器与数据集
void visual_pca(cv::Mat data,      // 可视化数据集，必须2两个特征
                cv::Mat &out,       // 输出的可视化图
                int w,          // 默认可视化区域
                int h){         // 默认可视化区域
    out = cv::Mat(h, w, CV_8UC3, cv::Scalar(255, 255, 255));
    // 计算缩放比例
    double scale_x, scale_y; 
    // 计算数据集的x，y最大最小值
    double max_x, min_x;
    double max_y, min_y;
    cv::minMaxIdx(data.col(0), &min_x, &max_x);
    cv::minMaxIdx(data.col(1), &min_y, &max_y);
    // 扩边，避免点绘制在边界上
    min_x -= 0.5;
    max_x += 0.5;
    min_y -= 0.5;
    max_y += 0.5;

    scale_x = w / (max_x - min_x);
    scale_y = h / (max_y - min_y);
    // 定义绘制的颜色
    cv::Vec3b r(0, 0, 255), g(0, 255, 0), b(255, 0, 0);
    // 绘制样本点
    for(int i = 0; i < data.rows; i++){
        double d_x = data.at<float>(i, 0);
        double d_y = data.at<float>(i, 1);
        // 先平移到0点再放大
        d_x -= min_x;
        d_y -= min_y;
        d_x *= scale_x;
        d_y *= scale_y; 
        // 取整绘制
        circle(out, cv::Point( (int) d_x,  (int) d_y ), 1, r, -1);
    }

    // 计算主成份
    cv::PCA pca(data, cv::Mat(), cv::PCA::DATA_AS_ROW);
    // 获取中心点
    double c_x = pca.mean.at<float>(0, 0);
    double c_y = pca.mean.at<float>(0, 1);
    c_x -= min_x;
    c_y -= min_y;
    c_x *= scale_x;
    c_y *= scale_y; 
    // 平移方法到（520，520）区域
    cv::Point cntr = cv::Point(
        static_cast<int>(c_x),
        static_cast<int>(c_y)
    );
    // 绘制中心点（均值点）
    circle(out, cntr, 5, cv::Scalar(255, 0, 255), 2);

    // 获取特征值与特征向量
    std::vector<cv::Point2d> eigen_vecs(2);
    std::vector<double> eigen_val(2);
    for (int i = 0; i < 2; i++){
        // 特征向量
        double e_x = pca.eigenvectors.at<float>(i, 0);
        double e_y = pca.eigenvectors.at<float>(i, 1);
        // 主成分向量
        eigen_vecs[i] = cv::Point2d(e_x,e_y);
        // 特征值(不放大改变)
        eigen_val[i] = pca.eigenvalues.at<float>(i);
    }
    
    // 计算两个主成分向量相对均值点的坐标(平移到均值点为中心):主成分需要变换后根据均值点平移
    cv::Point p1 = cntr + cv::Point(
        static_cast<int>(eigen_vecs[0].x * eigen_val[0] * scale_x), 
        static_cast<int>(eigen_vecs[0].y * eigen_val[0] * scale_y));
    cv::Point p2 = cntr - cv::Point(                            // 加减仅仅改变主成分向量相对均值点的方向
        static_cast<int>(eigen_vecs[1].x * eigen_val[1] * scale_x), 
        static_cast<int>(eigen_vecs[1].y * eigen_val[1] * scale_y));
    // 绘制主成分向量
    
    drawAxis(out, cntr, p1, cv::Scalar(0, 255, 0), 1.5);   // 放大1.5倍
    drawAxis(out, cntr, p2, cv::Scalar(255, 255, 0), 5.0); // 放大5倍
}

二.物体方向判别可视化

效果

物体方向判别可视化效果

代码

#include <opencv2/opencv.hpp>
#include <iostream>
#include <cmath>
// 函数定义
void drawAxis(cv::Mat&, cv::Point, cv::Point, cv::Scalar, const float);
// 使用PCA计算方向
double getOrientation(const std::vector<cv::Point> &, cv::Mat&);

int main(int argc, const char** argv) {
    // 读取图像
    cv::Mat src = cv::imread("../pca.jpg");
    // 使用灰度图像来处理，所以先转换为灰度图像
    cv::Mat gray;
    cv::cvtColor(src, gray, cv::COLOR_BGR2GRAY);
    // 把图像分成2值矩阵0-1
    cv::Mat bw;
    cv::threshold(gray, bw, 50, 255, cv::THRESH_BINARY | cv::THRESH_OTSU);

    // 按照轮廓，把图像分成多个数据集
    std::vector<std::vector<cv::Point> > contours;
    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;

        // 在原图像上绘制轮廓
        drawContours(src, contours, static_cast<int>(i), cv::Scalar(0, 0, 255), 2);
        // 绘制每个数据集的主成分到源图像。
        getOrientation(contours[i], src);

        cv::imwrite("../pca_out.png", src);

    }

    return 0;
}

void drawAxis(cv::Mat& img, cv::Point p, cv::Point q, cv::Scalar colour, const float scale = 0.2){
    /*
        cv::Mat& img        ：需要绘制的图像
        cv::Point p         ：绘制坐标轴线段的两个点之一：坐标轴的起点
        cv::Point q         ：绘制坐标轴线段的两个点之二：坐标轴的终点
        cv::Scalar colour   ：绘制的颜色
        const float scale   ：坐标轴的放大因子
     */

    // 计算两个点与原点形成的向量的夹角（用于延长线段的计算）
    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延长,并绘制（p不动，按照原来pq的联线方向，增加线段长度）
    q.x = (int) (p.x - scale * hypotenuse * cos(angle));
    q.y = (int) (p.y - scale * hypotenuse * sin(angle));
    line(img, p, q, colour, 1, cv::LINE_AA);

    // 绘制坐标轴的箭头（以终点q为起点的两个小线段，与原来线段成45°夹角），箭头线段9像素
    int len = 9;
    // 箭头：第一个线段(逆时针旋转45°)
    p.x = (int) (q.x + len * cos(angle + CV_PI / 4));
    p.y = (int) (q.y + len * sin(angle + CV_PI / 4));
    line(img, p, q, colour, 1, cv::LINE_AA);
    // 箭头：第二个线段(顺时针旋转45°)
    p.x = (int) (q.x + len * cos(angle - CV_PI / 4));
    p.y = (int) (q.y + len * sin(angle - CV_PI / 4));
    line(img, p, q, colour, 1, cv::LINE_AA);
}

// 计算图像中PCA降维向量
double getOrientation(const std::vector<cv::Point> &pts, cv::Mat &img){
    /*
        std::vector<cv::Point> &pts     ：数据集（从图像处理中获取）
        cv::Mat &img                    ：显示数据集的图像
     */

    // 先把std::vector<cv::Point> &pts 转换为cv::Mat，这样便于操作
    int sz = static_cast<int>(pts.size());
    cv::Mat data_pts(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降维
    cv::PCA pca(data_pts, cv::Mat(), cv::PCA::DATA_AS_ROW);  // 第二个参数为空，表示均值从第一个参数指定的数据计算。
    // 均值就是数据集的中心点
    cv::Point cntr = cv::Point(
        static_cast<int>(pca.mean.at<double>(0, 0)),    // 均值点第一个坐标
        static_cast<int>(pca.mean.at<double>(0, 1)));   // 均值点第二个坐标

    // 特征向量与特征值-因为是二维，所以是两个特征向量与特征值（特征值有可能是0）
    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.eigenvectors.at<double>(i, 0),
            pca.eigenvectors.at<double>(i, 1)
        );
        // 特征值
        eigen_val[i] = pca.eigenvalues.at<double>(i);
    }
    // 开始绘制
    // 1. 绘制中心点
    circle(img, cntr, 3, cv::Scalar(255, 0, 255), 2);
    // 2. 绘制第一个特征向量（主成分）
    // 计算绘制点
    cv::Point p1 = cntr + 0.02 * cv::Point(         // 像素太大，所以缩小到0.2
        static_cast<int>(eigen_vecs[0].x * eigen_val[0]), 
        static_cast<int>(eigen_vecs[0].y * eigen_val[0])
    );
    // 绘制
    drawAxis(img, cntr, p1, cv::Scalar(0, 255, 0), 1);   // 放大1倍(因为第一个主成分本身就很大，不再放大)

    // 3. 绘制第二个特征向量（主成分）
    cv::Point p2 = cntr - 0.02 * cv::Point(
        static_cast<int>(eigen_vecs[1].x * eigen_val[1]), 
        static_cast<int>(eigen_vecs[1].y * eigen_val[1])
    );
    drawAxis(img, cntr, p2, cv::Scalar(255, 255, 0), 5);  // 放大5倍（第二个特征向量太小，随意放大5倍）
    // 返回方向：第一个特征向量的方向
    double angle = atan2(eigen_vecs[0].y, eigen_vecs[0].x); 
    return angle;
}

cmake编译脚本

# 限制cmake的版本
cmake_minimum_required(VERSION 3.15)

# 设置C++标准的版本，OpenCV需要C++11支持
set(CMAKE_CXX_STANDARD 11)
set(CMAKE_CXX_STANDARD_REQUIRED True)

# 项目名
project(PCA)

# 加载OpenCV提供的cmake变量OpenCV_LIBS等，可以阅读OpenCV的cmake文件查阅其中的cmake变量清单
# set(CMAKE_PREFIX_PATH /usr/local)   
find_package(OpenCV REQUIRED)

# 设置编译链接的执行文件
add_executable(pca pca_img.cpp)
# 编译需要的动态库 或者 framework（Mac OS）
target_link_libraries(pca ${OpenCV_LIBS})

# #################第二个程序###########
add_executable(iris pca_iris.cpp)
target_link_libraries(iris ${OpenCV_LIBS})

---