PCANet的C++代码——详细注释版

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;
}