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