FastICA opencv 实现版

FastICA   opencv 实现版

Matlab版可参见：利用ICA进行图像加密

#include "stdafx.h"
#include <opencv.hpp>
#include <iostream>
using namespace cv;
using namespace std;


void fastIca(Mat S, int len1, int  len2)
{
	S = S/255;
	Mat MixedS = S.clone();
	Mat MixedS_bak = S.clone();

	int width = S.size().width;
	int heigh = S.size().height;

	Mat MixedS_mean(heigh,1,CV_64FC1);
	for (int i = 0; i < heigh; i++)
	{
		Mat MixS = MixedS.rowRange(i,i+1);
		double m = mean(MixS)[0];
		MixS = MixS - m;
	}

白化
	Mat MixedS_(MixedS.t());
	Mat covMat;  
	Mat meanMat;  
	calcCovarMatrix(MixedS_,covMat,meanMat,CV_COVAR_NORMAL|CV_COVAR_ROWS); //计算协方差
	Mat E,D;
	eigen(covMat, D, E);//计算特征值
	Mat eye= Mat::eye( E.size().height, E.size().height,CV_64FC1);
	for (int i = 0; i < E.size().height; i++)
	{
		Mat eyei = eye.rowRange(i,i+1);
		eyei = eyei*D.at<double>(i,0);
	}//单独特征值转换成矩阵
	pow(eye,0.5,eye);// eye = sqrt(D)
	Mat Q = eye.inv()*E.t();// Q = inv(sqrt(D))*E'
	Q.convertTo(Q,CV_64FC1);
	Mat MixedS_white = Q*MixedS;

FastICA//

	
	Mat X = MixedS_white.clone();
	int VariableNum = X.size().height;
	int SampleNum = X.size().width;
	int numofIC = VariableNum;
	Mat B = Mat::zeros(numofIC, VariableNum, CV_64FC1);
	Mat b(numofIC, 1, CV_64FC1);
	for (int i = 0; i < numofIC; i++)
	{
		int ii = 1;
		int maxIterationsNum = 100;
		int IterationsNum  = 0;
		randn(b,Scalar::all(0.2), Scalar::all(0.8)); //随机矩阵
		b = b - 0.5;
		b = b/norm(b);
		while (ii<maxIterationsNum)
		{
			if (ii == maxIterationsNum-1)
			{
				cout<<"第"<<ii<<"分量在"<<maxIterationsNum<<"次迭代内并不收敛"<<endl;
			}

			Mat b_Old = b;
			double a2 = 0.9;
			double u = 0.9;
			Mat t = X.t()*b;// t
			Mat tt(t.size(),CV_64FC1);
			pow(t,2,tt);// t.^2
			Mat tt_a(t.size(),CV_64FC1); 
			tt_a = -tt/2;// -t.^2/2
			exp(tt_a,tt_a); //tt_a = exp(-t.^2/2)
			Mat g = t.mul(tt_a);//g = t.*exp(-t.^2/2)
			Mat dg = (1-a2*tt).mul(tt_a);//dg = (1-a2*t.^2).*exp(-t.^2/2)
			Mat bb = (1-u)*t.t()*g; 
			double bb_ = bb.at<double>(0,0);
			Mat b_ = bb_*b;
			b_ = b_ +(u*X*g)/SampleNum; 
			b = b_-(mean(dg)[0]*b);   // b_ =((1-u)*t'*g*b+u*X*g)/SampleNum-mean(dg)*b;
			b = b-(B*B.t())*b;
			b = b/norm(b);
			Mat dst = b.t()*b_Old;
			double dstV = dst.at<double>(0,0);
			double flag = abs(abs(dstV)-1) ;
			if (flag  < 0.0000000000001   )
			{
					for (int j = 0; j < numofIC; j++)
					{
						B.at<double>(j,i) = b.at<double>(j,0);
						cout<<b<<endl;
						cout<<B<<endl;
					}
					break;
			}
			++ii;
		}
	}
	Mat ICAedS = B.t()*Q*MixedS_bak;
	double maxO , minO;
	minMaxIdx(abs(ICAedS), &minO, &maxO);
	ICAedS = abs(ICAedS)/maxO;
	for (int i = 0; i < heigh; i++)
	{
		Mat out1 = ICAedS.rowRange(i,i+1);
		Mat out2 = out1.reshape(1,len2);
		double maxV, minV;
		minMaxIdx(out2,&minV,&maxV);
		out2 = out2/maxV;
		out2 = out2 *255;
		out2.convertTo(out2,CV_8UC1);
		stringstream ss;
		ss <<i;
		imshow(ss.str(),out2);
		waitKey();
	}

}
void recoverB(Mat &img , int wid, int hei)
{
	int width = img.size().width;
	int height = img.size().height;
		for (int i = 0; i < height; i++)
	{
		Mat out1 = img.rowRange(i,i+1);
		Mat out2 = out1.reshape(1,hei);
		double maxV, minV;
		minMaxIdx(out2,&minV,&maxV);
		out2 = out2/maxV;
		out2 = out2 *255;
		out2.convertTo(out2,CV_8UC1);
		String u() ;
		stringstream ss;
		ss <<i;
		imshow(ss.str(),out2);
		waitKey();
	}

}

Mat mergeRows(Mat A, Mat B)
{
    int totalRows = A.rows + B.rows;

    Mat mergedDescriptors(totalRows, A.cols, A.type());
    Mat submat = mergedDescriptors.rowRange(0, A.rows);
    A.copyTo(submat);
    submat = mergedDescriptors.rowRange(A.rows, totalRows);
    B.copyTo(submat);
    return mergedDescriptors;
}
int _tmain(int argc, _TCHAR* argv[])
{

	Mat src1 = imread("src1.png",0);//读取两幅线性混叠的图像或者语音信号
	Mat src2 = imread("src2.jpg",0);//（如果是语音信号，reshape部分需要更改）
	double maxV1,minV1; 
	minMaxIdx(src1,&minV1,&maxV1);
	double maxV2, minV2;
	minMaxIdx(src2,&minV2,&maxV2);

	int len1 = src1.size().area();
	int len2 = src2.size().area();

	Mat dst1 = src1.reshape(1,1);//第一幅图转化成1*N维向量
	Mat dst2 = src2.reshape(1,1);//第二幅图转化成1*N维向量
	
	Mat dst_all = mergeRows(dst1,dst2);//将第一幅和第二幅图合并
	Mat weight(dst_all.size().height,dst_all.size().height,CV_64FC1);
	randu(weight,0,0.5);//随机权重
	dst_all.convertTo(dst_all,CV_64FC1);//
	dst_all  = dst_all/255;
	dst_all = weight * dst_all;//乘以随机权重后的矩阵

	recoverB(dst_all,src1.size().width, src1.size().height);
	//用于查看混淆后的样子

	fastIca( dst_all , src1.size().width, src1.size().height);
	//FastICA 算法实现


	return 0;
}