線性代數的直覺理解（5）

抱歉，文章許久沒有更新。接下來，這裡將要再次活躍起來！

PCA降維在數字圖像處理中的應用

上一篇文章的後半部分，我們提到了PCA的數學原理。這次來點兒實戰：結合OpenCV在

其核心功能模塊（core functionality)當中的PCA類，直觀看看數字圖像處理當中PCA的應用。

近幾年流行的代碼託管網站github上面有我不定期更新的代碼哦~歡迎去看熱鬧，更歡迎點那個小星星！

Detailed Description

Principal Component Analysis.

The class is used to calculate a special basis for a set of vectors. The basis will consist of eigenvectors of the covariance matrix calculated from the input set of vectors. The class PCA can also transform vectors to/from the new coordinate space defined by the basis. Usually, in this new coordinate system, each vector from the original set (and any linear combination of such vectors) can be quite accurately approximated by taking its first few components, corresponding to the eigenvectors of the largest eigenvalues of the covariance matrix. Geometrically it means that you calculate a projection of the vector to a subspace formed by a few eigenvectors corresponding to the dominant eigenvalues of the covariance matrix. And usually such a projection is very close to the original vector. So, you can represent the original vector from a high-dimensional space with a much shorter vector consisting of the projected vectors coordinates in the subspace. Such a transformation is also known as Karhunen-Loeve Transform, or KLT. See http://en.wikipedia.org/wiki/Principal_component_analysis

簡單地概括一下，就是這個類能夠計算出一組向量空間下最優的一組基向量，有了這組基向量，就很容易地調用裡面的公有方法完成向該基的投影求出新的坐標值，或者是能夠將新基下的坐標反向投影求出其原向量。我們首先看看這個類的屬性，屬性決定了這個類到底有多強大！

Public Attributes

Mat eigenvalues eigenvalues of the covariation matrix More...

Mat eigenvectors eigenvectors of the covariation matrix More...

Mat mean mean value subtracted before the projection and added after the back projection More...

有特徵值，特徵向量以及在求協方差矩陣之前一步必須中心化的那個矩陣。有了這三個關鍵點，我想這個類的成員方法應該差不多清楚了吧：

PCA () default constructor More...

PCA (InputArray data, InputArray mean, int flags, int maxComponents=0) PCA (InputArray data, InputArray mean, int flags, double retainedVariance) Mat backProject (InputArray vec) const Reconstructs vectors from their PC projections. More...

void backProject (InputArray vec, OutputArray result) const PCA & operator() (InputArray data, InputArray mean, int flags, int maxComponents=0) performs PCA More...

PCA & operator() (InputArray data, InputArray mean, int flags, double retainedVariance) Mat project (InputArray vec) const Projects vector(s) to the principal component subspace. More...

void project (InputArray vec, OutputArray result) const void read (const FileNode &fn) load PCA objects More...

void write (FileStorage &fs) const write PCA objects More...

構造函數，一般我們最好使用帶參數的構造函數對這個類進行實例化。上面兩個帶參數的構造函數，區別就是最後一個參數：一個是保留的最大維數值，默認是全部保留，另一個的描述是Percentage of variance that PCA should retain，即經過PCA變換之後，保留前百分之多少的特徵向量作為基底，這種初始化方法能夠做到歸一化處理，避免由於輸入數據的維度不同產生不同的特徵向量最大值。其實不難發現，當構造函數完成之後，那組最優的基底已經計算出來。剩下的方法例如

projectnbackProjectn

project是基於基底的投影或者backProject是給定係數求基底的線性組合了。

結合官方文檔，下面是代碼部分：

程序主要分成如下幾個步驟：

1. 讀取圖像序列
2. PCA準備階段：「拉直」圖像，構成扁平矩陣
3. 執行PCA
4. 得到圖像向量協方差矩陣的特徵向量之後，投影產生point，圖像「復原」
5. 使用trackbar動態調整保留特徵向量數量比例

下面上代碼：

#include <iostream>n#include <fstream>n#include <sstream>n#include <opencv2/core.hpp>n#include "opencv2/imgcodecs.hpp"n#include <opencv2/highgui.hpp>n#include "cvui.h"nnusing namespace cv;nusing namespace std;nn#define WINDOW_NAME "PCA"nint main(int argc, char** argv)n{n cv::CommandLineParser parser(argc, argv, keys);n string file = parser.get<string>("filename");n vector<Mat> imgs;n tryn {n readImgList(file,imgs);n }n catch(cv::Exception& e)n {n cerr<<"Error opening file ""<<file<<"". Reason:"<<e.msg<<endl;n exit(1);n }n //this demo needs at least 2 imagesn if(imgs.size()<=1)n {n string error_message = "This demo needs at least 2 images to work! Please add more image to your data set!";n CV_Error(Error::StsError,error_message);n }n //reshape and stack images into a row matrixn Mat data = formatImagesForPCA(imgs);n //perform pcan PCA pca(data,cv::Mat(),PCA::DATA_AS_ROW,0.95);n //demostration of the effect of retainedVaraince on the first imagen Mat point = pca.project(data.row(0));//project into the eigenspace,thus the image becomes a point with each elementn //representing the linear combination coefficientsn Mat reconstruction = pca.backProject(point);//recreate the image from the "point"n reconstruction = reconstruction.reshape(imgs[0].channels(),imgs[0].rows);n reconstruction = toGrayscale(reconstruction);nn Mat frame = Mat(80,300,CV_8UC3);n double trackBarValue = 100.0;n bool checked1 = false;n bool checked2 = true;n namedWindow(WINDOW_NAME);n cvui::init(WINDOW_NAME);n namedWindow("Face",WINDOW_AUTOSIZE);n for(;;)n {n //fill the frame with a nice color-RGBn frame = Scalar(40,50,40);n //show some textn cvui::text(frame,10,5,"This is PCA demostration");n //trackBar componentn cvui::trackbar(frame,50,30,200,&trackBarValue,0.,100.);nn PCA pca_realTime(data,cv::Mat(),PCA::DATA_AS_ROW,trackBarValue/100.0);n //demostration of the effect of retainedVaraince on the first imagen Mat point_realTime = pca_realTime.project(data.row(0));//project into the eigenspace,thus the image becomes a point with each elementn //representing the linear combination coefficientsn Mat reconstruction_realTime = pca_realTime.backProject(point_realTime);//recreate the image from the "point"n reconstruction_realTime = reconstruction_realTime.reshape(imgs[0].channels(),imgs[0].rows);n reconstruction_realTime = toGrayscale(reconstruction_realTime);n //update must be called after all componentsn //it handles mouse click and others behind the scene with magic!n cvui::update();n imshow(WINDOW_NAME,frame);n imshow("Face",reconstruction_realTime);n if(waitKey(20)==27)//press escp to exit the processn {n break;n }n }nn return 0;n}n

1.讀取數據源

進入主函數之後，首先通過CommandLineParser對象方便地讀取用戶指定文件夾下的圖像文件。由於後續的人臉識別演算法也要用到data文件夾下的圖像，但是本次程序只用到了每個文件夾下的1.pgm文件，其他的都是同一個人不同角度獲取的圖像，暫不使用。我的工程結構是這樣的：

/rootn main.cppn imgList.txtn /buildn Pca Unix Excutablen /datan /s1/1.pgmn /s1/2.pgmn …… n /s2n ……n

imgList.txt文件，說明了readImgList方法的文件路徑，下面分別是路徑以及方法定義：

../data/s1/1.pgmn../data/s2/1.pgmn../data/s3/1.pgmn../data/s4/1.pgmn../data/s5/1.pgmn../data/s6/1.pgmn../data/s7/1.pgmn../data/s8/1.pgmn../data/s9/1.pgmn../data/s10/1.pgmn../data/s11/1.pgmn../data/s12/1.pgmn../data/s13/1.pgmn../data/s14/1.pgmn../data/s15/1.pgmn

readImgList：

static void readImgList(const string& filename,vector<Mat>& images)n{n std::ifstream file(filename.c_str(), ifstream::in);n if (!file) {n string error_message = "No valid input file was given, please check the given filename.";n CV_Error(Error::StsBadArg, error_message);n }n string line;n while (getline(file, line)) {n images.push_back(imread(line, 0));n }n}n

2.拉直

圖像對象存儲在了imgs容器中，接下來就開始第二步驟，圖像「拉直」。

//each image M*N,then the output Mat is numberof(images)*(M*N) flat matrixnstatic Mat formatImagesForPCA(const vector<Mat>& data)n{n Mat dst(static_cast<int>(data.size()),data[0].rows*data[0].cols,CV_32F);n for(unsigned int i=0;i<data.size();i++)n {n Mat image_row = data[i].clone().reshape(1,1);n Mat row_i = dst.row(i);n image_row.convertTo(row_i,CV_32F);n }n return dst;n}n

執行上面的方法之後，扁平矩陣data構造完成：

Mat data = formatImagesForPCA(imgs);n

3.執行PCA

計算輸入數據協方差矩陣的特徵向量，保留前面95%。

//perform pcanPCA pca(data,cv::Mat(),PCA::DATA_AS_ROW,0.95);n

注意，這裡第二個參數，默認設置成為empty，PCA對象首先自動中心化，然後再計算協方差矩陣，求特徵向量。

mean optional mean value; if the matrix is empty (noArray()), the mean is computed from the data.n

4.投影和圖像「復原」

//demostration of the effect of retainedVaraince on the first imagenMat point = pca.project(data.row(0));//project into the eigenspace,thus the image becomes a point with each elementn//representing the linear combination coefficientsnMat reconstruction = pca.backProject(point);//recreate the image from the "point"nreconstruction = reconstruction.reshape(imgs[0].channels(),imgs[0].rows);nreconstruction = toGrayscale(reconstruction);n

point是將第一幅拉直後的圖像投影得到的新基下向量。例如，本例中，單張灰度圖解析度為112*92，那麼拉直之後，就是10304的列向量。調試運行程序，可以看到特徵向量按行排列，共15行，有15個特徵向量產生。

上面的調試結果是PCA構造函數最後一個參數設置100%得到的，當設置成95%時，計算產生了11個特徵向量值。經point重新復原的圖像，調用了toGrayscale方法，目的是將各個像素點灰度值比例限定在0~255之間，為顯示圖像做準備。

static Mat toGrayscale(InputArray _src)n{n Mat src = _src.getMat();n //only allow one channeln if(src.channels()!=1)n {n CV_Error(Error::StsBadArg,"Only matrices with one channel are supported!");n }n //create and return normolized imagen Mat dst;n cv::normalize(_src,dst,0,255,NORM_MINMAX,CV_8UC1);n return dst;n}n

5.使用trackbar動態調整特徵向量保留比例

這裡使用了一種簡易的UI界面交互小控制項，使用方法能夠在我的倉庫中找到。

完成上述過程，得到如下的界面：