PCANet的C++代码——详细注释版
2016-04-21 21:46
441 查看
PCANet的C++代码:here。
代码注释:
代码注释:
cv::Mat im2colstep(cv::Mat& InImg, vector<int>& blockSize, vector<int>&stepSize) { // block的行大小 * 列大小,定义OutBlock的列数。 int r_row = blockSize[ROW_DIM] * blockSize[COL_DIM]; // 图片的rows减去block的rows, 为了防止在滑动窗口的时候越界。 图片的发小要减去block的大小。 int row_diff = InImg.rows - blockSize[ROW_DIM]; int col_diff = InImg.cols - blockSize[COL_DIM]; // 减完之后要除以步长需要再加1。计算所有block的个数,作为OutBlock的行数。 int r_col = (row_diff / stepSize[ROW_DIM] + 1) * (col_diff / stepSize[COL_DIM] + 1); // 定义一个Mat 类型,用于保存提取的所有block。 cv::Mat OutBlocks(r_col, r_row, InImg.depth()); // 定义指向OutBlocks的行指针。 double *p_InImg, *p_OutImg; // 计数。当前第几个block,或者说OutBlocks的第几行。 int blocknum = 0; // block横向遍历图片。 for(int j=0; j<=col_diff; j+=stepSize[COL_DIM]) { // block纵向遍历图片。 for(int i=0; i<row_diff; i+=stepSize[ROW_DIM]) { // 初始化指针。 p_OutImg = OutBlocks.ptr<double>(blocknum); // 下面就是将每个block的值保存到p_OutImg中。 // 用来计数行数。 for(int m=0; m<blockSize[ROW_DIM]; m++) { // 将InImg的第(i+m)行赋值给p_InImg。 p_InImg = InImg.ptr<double>(i + m); // 用来计数列数。 for(int l=0; l<blockSize[COL_DIM]; l++) { // 将p_InImg数组中的j+l个值赋给p_OutImg。也就是block的m行l列。 // 注意在存放的时候是以列优先。 p_OutImg[blockSize[ROW_DIM] * l + m] = p_InImg[j+l]; } } blocknum++; } } return OutBlocks; }
// 得到图片的所有channels的所有block。 cv::Mat im2col_general(cv::Mat& InImg, vector<int>& blockSize, vector<int>& stepSize) { // block的大小,表示行和列。 & 遍历图片行列的步长。 assert(blockSize.size() == 2 && stepSize.size() == 2); // 图片的通道数。 int channels = InImg.channels(); // 定义一个向量,保存InImg图片的每个通道。 vector<cv::Mat> layers; // 如果InImg的通道数大于1,需要分割。 if(channels > 1) split(InImg, layers); else // 通道数为1,不需要再分割。 layers.push_back(InImg); // 得到InImg第一个通道图片的所有block。 cv::Mat AllBlocks = im2colstep(layers[0], blockSize, stepSize); int size = layers.size(); if(size > 1) { // 交换layers第一个元素和最后一个元素。目的:删除layers的第一个元素。 swap(layers[0], layers.back()); layers.pop_back(); for(int i=1; i<size; i++) { // 实现图像的拼接。 将InImg的每个通道拼接到一起。 hconcat(AllBlocks, im2colstep(layers[i], blockSize, stepSize)); } } // 矩阵的转置。 return AllBlocks.t(); }
// 求样本的特征向量。 cv::Mat PCA_FilterBank(vector<cv::Mat>& InImg, int PatchSize, int NumFileters) { int channels = InImg[0].channels(); int InImg_Size = InImg.size(); int* randIdx = getRandom(InImg_Size); // size用来定义协方差矩阵的大小。 int size = channels * PatchSize * PatchSize; int img_depth = InImg[0].depth(); cv::Mat Rx = cv::Mat::zeros(size, size, img_depth); // 设置block的行数和列数。 设置横向遍历和纵向遍历的步长。 vector<int> blockSize; vector<int> stepSize; for(int i=0; i<2; i++) { // 可见选取的block长宽相等,步长都为1。 blockSize.push_back(PatchSize); stepSize.push_back(1); } // 用于保存得到的block。 cv::Mat temp; // 存放样本的每个特征的均值。 cv::Mat mean; // 存放样本每维减去对应的均值。 cv::Mat temp2; cv::Mat temp3; // OpenMP指令。 int coreNum = omp_get_num_procs(); int cols = 0; # pragma omp parallel for default(none) num_threads(coreNum) private(temp, temp2, temp3, mean) shared(cols, Rx, InImg_Size, InImg, randIdx, blockSize, stepSize) for(int j=0; j<InImg_Size; j++) { temp = im2col_general(InImg[randIdx[j]], blockSize, stepSize); // PCA第一步,求特征均值。在这里是以列为特征的,下边为什么要以行为特征? cv::reduce(temp, mean, 0, CV_REDUCE_AVG); temp3.create(0,temp.cols, temp.type()); cols = temp.cols; // 所有样本减去均值。 for(int i=0; i<temp.rows; i++) { temp2 = (temp.row(i) - mean.row(0)); temp3.push_back(temp2.row(0)); } temp = temp3 * temp3.t(); # pragma omp critical Rx = Rx + temp; } Rx = Rx / (double)(InImg_Size * cols); cv::Mat eValuesMat; cv::Mat eVectorsMat; // 求特征值和特征向量。 eigen(Rx, eValuesMat,eVectorsMat); cv::Mat Filters(0, Rx.cols, Rx.depth()); // 选取NumFileters个特征向量作为行向量组成特征向量矩阵。 for(int i=0; i<NumFileters; i++) Filters.push_back(eVectorsMat.row(i)); return Filters; }
// 接下来,减去均值样本矩阵 * 特征向量就是PCA的输出。 PCA_Out_Result* PCA_output(vector<cv::Mat>& InImg, vector<int>& InImgIdx, int PatchSize, int NumFilters, cv::Mat& Filters, int threadnum){ PCA_Out_Result* result = new PCA_Out_Result; int img_length = InImg.size(); int mag = (PatchSize - 1) / 2; int channels = InImg[0].channels(); cv::Mat img; vector<int> blockSize; vector<int> stepSize; for(int i=0; i<2; i++){ blockSize.push_back(PatchSize); stepSize.push_back(1); } cv::Mat temp; cv::Mat mean; cv::Mat temp2; cv::Mat temp3; int coreNum = omp_get_num_procs();//»ñµÃŽŠÀíÆ÷žöÊý coreNum = coreNum > threadnum ? threadnum : coreNum; cv::Scalar s = cv::Scalar(0); # pragma omp parallel for default(none) num_threads(coreNum) private(img, temp, temp2, temp3, mean) shared(InImgIdx, s, blockSize, stepSize, mag, img_length, InImg, result, Filters, NumFilters) for(int i=0; i<img_length; i++){ //对样本进行边缘补零操作,以保证映射结果与原图像的尺寸相同。 cv::copyMakeBorder(InImg[i], img, mag, mag, mag, mag, cv::BORDER_CONSTANT, s); temp = im2col_general(img, blockSize, stepSize); cv::reduce(temp, mean, 0, CV_REDUCE_AVG); temp3.create(0, temp.cols, temp.type()); for(int i=0;i<temp.rows;i++){ temp2 = (temp.row(i) - mean.row(0)); temp3.push_back(temp2.row(0)); } #pragma omp critical { // 保存图片的索引号。 result->OutImgIdx.push_back(InImgIdx[i]); for(int j=0; j<NumFilters; j++) { // 减去均值样本矩阵 * 特征向量 temp = Filters.row(j) * temp3; temp = temp.reshape(0, InImg[i].cols); result->OutImg.push_back(temp.t()); } } } /* int size = InImgIdx.size(); for(int i=0; i<size; i++) for(int j=0; j<NumFilters; j++) result->OutImgIdx.push_back(InImgIdx[i]);*/ return result; }
PCA_Train_Result* PCANet_train(vector<cv::Mat>& InImg, PCANet* PcaNet, bool is_extract_feature){ assert(PcaNet->NumFilters.size() == PcaNet->NumStages); PCA_Train_Result* train_result = new PCA_Train_Result; PCA_Out_Result* out_result = new PCA_Out_Result; PCA_Out_Result* temp; out_result->OutImg = InImg; int img_length = InImg.size(); for(int i=0; i<img_length; i++) out_result->OutImgIdx.push_back(i); int64 e1 = cv::getTickCount(); int64 eo1, eo2, eb1, eb2; for(int s=0; s<PcaNet->NumStages; s++){ eb1 = cv::getTickCount(); cout << " Computing PCA filter bank and its outputs at stage " << s << "..." << endl; train_result->Filters.push_back(PCA_FilterBank(out_result->OutImg, PcaNet->PatchSize, PcaNet->NumFilters[s])); eb2 = cv::getTickCount(); cout <<" stage"<<s<<" PCA_FilterBank time: "<<(eb2 - eb1)/ cv::getTickFrequency()<<endl; eo1 = cv::getTickCount(); if(s != PcaNet->NumStages - 1){ temp = PCA_output(out_result->OutImg, out_result->OutImgIdx, PcaNet->PatchSize, PcaNet->NumFilters[s], train_result->Filters[s], omp_get_num_procs()); delete out_result; out_result = temp; } eo2 = cv::getTickCount(); cout <<" stage"<<s<<" output time: "<<(eo2 - eo1)/ cv::getTickFrequency()<<endl; } int64 e2 = cv::getTickCount(); double time = (e2 - e1)/ cv::getTickFrequency(); cout <<"\n totle FilterBank time: "<<time<<endl; InImg.clear(); vector<cv::Mat>().swap(InImg); vector<cv::Mat> tempF; int end = PcaNet->NumStages - 1; int outIdx_length = out_result->OutImgIdx.size(); if(is_extract_feature){ vector<cv::Mat>::const_iterator first = out_result->OutImg.begin(); vector<cv::Mat>::const_iterator last = out_result->OutImg.begin(); vector<cv::Mat> features; Hashing_Result* hashing_r; int coreNum = omp_get_num_procs();//»ñµÃŽŠÀíÆ÷žöÊý e1 = cv::getTickCount(); # pragma omp parallel for default(none) num_threads(coreNum) private(temp, hashing_r) shared(features, out_result, PcaNet, first, last, outIdx_length, img_length, train_result, end) for(int i=0; i<img_length; i++){ vector<cv::Mat> subInImg(first + i * PcaNet->NumFilters[end], last + (i + 1) * PcaNet->NumFilters[end]); vector<int> subIdx; /*for(int j=0; j<outIdx_length; j++){ if(out_result->OutImgIdx[j] == i) subIdx.push_back(1); else subIdx.push_back(0); }*/ for(int j=0; j< PcaNet->NumFilters[end]; j++) subIdx.push_back(j); temp = PCA_output(subInImg, subIdx, PcaNet->PatchSize, PcaNet->NumFilters[end], train_result->Filters[end], 2); hashing_r = HashingHist(PcaNet, temp->OutImgIdx, temp->OutImg); #pragma omp critical { features.push_back(hashing_r->Features); train_result->feature_idx.push_back(out_result->OutImgIdx[i]); } delete hashing_r; delete temp; subIdx.clear(); vector<int>().swap(subIdx); } e2 = cv::getTickCount(); time = (e2 - e1)/ cv::getTickFrequency(); cout <<"\n hasing time: "<<time<<endl; //out_result->OutImg.clear(); //vector<cv::Mat>().swap(out_result->OutImg); delete out_result; int size = features.size(); if(size > 0){ train_result->Features.create(0, features[0].cols, features[0].type()); for(int i=0 ;i<size; i++){ train_result->Features.push_back(features[i]); } /* train_result->Features = features[0]; for(int i=1 ;i<size; i++){ vconcat(train_result->Features, features[i], train_result->Features); }*/ } features.clear(); vector<cv::Mat>().swap(features); } //if(temp != NULL) // delete temp; return train_result; }
double round(double r) { return (r > 0.0) ? floor(r + 0.5) : ceil(r - 0.5); }
Hashing_Result* HashingHist(PCANet* PcaNet, vector<int>& ImgIdx, vector<cv::Mat>& Imgs) { Hashing_Result* ha_result = new Hashing_Result; int length = Imgs.size(); // 得到最后层的PCA的特征向量数。也就是第二阶段产生的图片数。 int NumFilters = PcaNet->NumFilters[PcaNet->NumStages - 1]; // 得到第一层PCA的特征向量数。也就是第一阶段产生的图片数。 int NumImgin0 = length / NumFilters; cv::Mat T; int row = Imgs[0].rows; int col = Imgs[0].cols; int depth = Imgs[0].depth(); // 这部分就是进行二值化哈希编码时候的权重。 vector<double> map_weights; cv::Mat temp; for(int i=NumFilters - 1; i>=0; i--) { // 权重设置为 2 的 i 次幂。 map_weights.push_back(pow(2.0,(double)i)); } // 保存 滑动窗口的横向步长和纵向步长。 vector<int> Ro_BlockSize; // 可以说成是滑动窗口的覆盖率吧。 double rate = 1 - PcaNet->BlkOverLapRatio; for(int i=0; i<PcaNet->HistBlockSize.size(); i++) { // HistBlockSize是滑动窗口的大小,包括宽和高。 Ro_BlockSize.push_back(round(PcaNet->HistBlockSize[i] * rate)); } cv::Mat BHist; int ImgIdx_length = ImgIdx.size(); int* new_idx = new int[ImgIdx_length]; for(int i=0; i<ImgIdx_length; i++) { new_idx[ImgIdx[i]] = i; } // PCA的输出一共有64个mat。 for(int i=0; i<NumImgin0; i++) { T = cv::Mat::zeros(row,col,depth); for(int j=0; j<NumFilters; j++) { // Heaviside进行二值化处理。 temp = Heaviside(Imgs[NumFilters * new_idx[i] + j]); // 乘以权重,每个像素值都被编码成为(0,255)之间的整数。 temp = temp * map_weights[j]; // 之后8个mat进行累加。 T = T + temp; } // 第一层的每个输出矩阵,将其分为B块,计算统计每个块的直方图信息,然后在将各个块的直方图特征进行级联,最终得到块扩展直方图特征. temp = im2col_general(T, PcaNet->HistBlockSize, Ro_BlockSize); temp = Hist(temp, (int)(pow(2.0, NumFilters)) - 1); temp = bsxfun_times(temp, NumFilters); if(i == 0) BHist = temp; else hconcat(BHist,temp,BHist); } int rows = BHist.rows; int cols = BHist.cols; ha_result->Features.create(1, rows * cols, BHist.type()); double *p_Fe = ha_result->Features.ptr<double>(0); double *p_Hi; for(int i=0; i<rows; i++) { p_Hi = BHist.ptr<double>(i); for(int j=0; j<cols; j++) { p_Fe[j * rows + i] = p_Hi[j]; } } return ha_result; }
// 进行二值化处理。 cv::Mat Heaviside(cv::Mat& X) { int row = X.rows; int col = X.cols; int depth = X.depth(); cv::Mat H(row, col, depth); double *p_X, *p_H; ////# pragma omp parallel for default(none) num_threads(4) private(p_X, p_H) shared(X, H, row, col) // 对于X的每一行中的每一列的元素,如果>0,则变为1;否则为0。 for(int i=0; i<row; i++) { p_X = X.ptr<double>(i); p_H = H.ptr<double>(i); for(int j=0; j<col; j++) { if(p_X[j]>0) p_H[j] = 1; else p_H[j] = 0; } } return H; }
// 获得直方图 cv::Mat Hist(cv::Mat& mat, int Range) { cv::Mat temp = mat.t(); int row = temp.rows; int col = temp.cols; int depth = temp.depth(); cv::Mat Hist = cv::Mat::zeros(row, Range + 1,depth); double *p_M, *p_H; ////# pragma omp parallel for default(none) num_threads(4) private(p_M, p_H) shared(temp, Hist, row, col) for(int i=0; i<row; i++) { p_M = temp.ptr<double>(i); p_H = Hist.ptr<double>(i); for(int j=0; j<col; j++) { p_H[(int)p_M[j]] += 1; } } temp = Hist.t(); return temp; }
// 直方图特征。 cv::Mat bsxfun_times(cv::Mat& BHist, int NumFilters) { double *p_BHist; int row = BHist.rows; int col = BHist.cols; double* sum = new double[col]; // 初始化sum。 for(int i=0; i<col; i++) { sum[i] = 0; } // 计算所有行元素的相加和。 for(int i=0; i<row; i++) { p_BHist = BHist.ptr<double>(i); for(int j=0; j<col; j++) { sum[j] += p_BHist[j]; } } double p = pow(2.0, NumFilters); ////# pragma omp parallel for default(none) num_threads(4) shared(col, sum, p) for(int i=0; i<col; i++) { sum[i] = p/sum[i]; } ////# pragma omp parallel for default(none) num_threads(4) private(p_BHist) shared(col, row, sum, BHist) // 对于每行元素乘以sum。 for(int i=0; i<row; i++) { p_BHist = BHist.ptr<double>(i); for(int j=0; j<col; j++) { p_BHist[j] *= sum[j]; } } return BHist; }
相关文章推荐
- c++第4次作业
- C++流操作算子
- 消失的字符串——c语言函数中的数据存储方式以及字符串实现
- 关于C语言的问卷调查
- Android NDK开发C语言部分的单步跟踪
- 第四周作业——C语言自评
- C语言
- C++ Greedy Snake的OOP实现 贪食蛇 <list> STL初次学习
- 两人合作审阅C++装饰模式
- 二叉排序树的C++实现,包括难点删除
- C语言贪食蛇
- 2016年4月21日 21:18:25 我的第一篇博客~
- effective C++ 读书笔记 条款08
- c和c++栈
- c++map基本操作
- 如何使用VC++6.0发布程序(即release版本程序)
- 程序设计篇(1):学生经验值管理系统(单链表实现)
- 山东省第五届ACM大学生程序设计竞赛-Weighted Median(模拟)
- 循环队列的综合操作(出入队、获取队长度)C语言实现
- char* p与char p[]