function [ dst ] = compute( I1,I2 ,src,ii )
%UNTITLED3 Summary of this function goes here
% Detailed explanation goes here
maxV1 = double(max(max(I1)))/255;
minV1 = double(min(min(I1)))/255;
maxV2 = double(max(max(I2)))/255;
minV2 = double(min(min(I2)))/255;
M = size(I1,1)*size(I1,2);
len1 = size(I1,1);
len2 = size(I1,2);
src1 = double(reshape(I1,1,M))/255;
src2 = double(reshape(I2,1,M))/255;
S=[src1;src2];
MixedS=S;
MixedS_bak=S;
%%%%%%%%%%%%%%%%%%%%%%%%%% 标准化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_mean=zeros(2,1);
for i=1:2
MixedS_mean(i)=mean(MixedS(i,:));
end % 计算MixedS的均值
for i=1:2
for j=1:size(MixedS,2)
MixedS(i,j)=MixedS(i,j)-MixedS_mean(i);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%% 白化 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MixedS_cov=cov(MixedS'); % cov为求协方差的函数
[E,D]=eig(MixedS_cov); % 对信号矩阵的协方差函数进行特征值分解
Q=inv(sqrt(D))*(E)'; % Q为白化矩阵
MixedS_white=Q*MixedS; % MixedS_white为白化后的信号矩阵
IsI=cov(MixedS_white'); % IsI应为单位阵
%%%%%%%%%%%%%%%%%%%%%%%% FASTICA算法 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X=MixedS_white; % 以下算法将对X进行操作
[VariableNum,SampleNum]=size(X);
numofIC=VariableNum; % 在此应用中，独立元个数等于变量个数
B=zeros(numofIC,VariableNum); % 初始化列向量w的寄存矩阵,B=[b1 b2 ... bd]
for r=1:numofIC
i=1;maxIterationsNum=100; % 设置最大迭代次数（即对于每个独立分量而言迭代均不超过此次数）
IterationsNum=0;
b=rand(numofIC,1)-.5; % 随机设置b初值
b=b/norm(b); % 对b标准化 norm(b):向量元素平方和开根号
while i<=maxIterationsNum+1
if i == maxIterationsNum % 循环结束处理
fprintf('\n第%d分量在%d次迭代内并不收敛。', r,maxIterationsNum);
break;
end
bOld=b;
a2=1;
u=1;
t=X'*b;
g=t.*exp(-a2*t.^2/2);
dg=(1-a2*t.^2).*exp(-a2*t.^2/2);
b=((1-u)*t'*g*b+u*X*g)/SampleNum-mean(dg)*b;
% 核心公式
b=b-B*B'*b; % 对b正交化
b=b/norm(b);
if abs(abs(b'*bOld)-1)<1e-9 % 如果收敛，则
B(:,r)=b; % 保存所得向量b
break;
end
i=i+1;
end
% B(:,r)=b; % 保存所得向量b
end
%%%%%%%%%%%%%%%%%%%%%%%%%% ICA计算的数据复原并构图 %%%%%%%%%%%%%%%%%%%%%%%%%
ICAedS=B'*Q*MixedS_bak; % 计算ICA后的矩阵
% 将混合矩阵重新排列并输出
ICAedS=abs(ICAedS)/max(max(abs(ICAedS)));
minv1 = min(ICAedS(1,:));
maxv1 = max(ICAedS(1,:));
minv2 = min(ICAedS(2,:));
maxv2 = max(ICAedS(2,:));
%% 错误提前判断
if minv1<minV1
dst1 = imadjust(reshape(ICAedS(1,:),len1,len2),[minv1,maxv1],[minV1,maxV1],1);
else
dst1 = imadjust(reshape(ICAedS(1,:),len1,len2));
end
if minv2<minV2
dst2 = imadjust(reshape(ICAedS(2,:),len1,len2),[minv2,maxv2],[minV2,maxV2],1);
else
dst2 = imadjust(reshape(ICAedS(2,:),len1,len2));
end
%% 得到分离图像
if sum(sum(dst2))> sum(sum(dst1))
dst1(dst1<0.12)=0;
dst = dst1;
else
dst2(dst2<0.12)=0;
dst = dst2;
end
[h w] = size(dst);
%% 获取要处理的区域
getRoi = dst(h*0.2:h*0.4,:);
dd = getRoi>0.2;%% 消除过暗的区域
[hh,ww] = size(getRoi);
imshow(imresize(src,0.5));
hold on
%% 画图
s = regionprops(dd);
flag = 0;
for r = 1:length(s)
if s(r).Area<20 || s(r).BoundingBox(2)>hh*0.8 ||s(r).BoundingBox(2)<hh*0.05|| s(r).BoundingBox(1)>ww*0.85 ||s(r).BoundingBox(1)<(ww*0.15)
continue
end
flag = 1;
% plot(s(r).Centroid(:,1), s(r).Centroid(:,2), 'r*') %%%中心点
%% 映射到原图
% uu = s(r).BoundingBox+[0,h*0.2,0,0];%%%映射到原图
% uu = uu*2;
% rectangle('Position',uu,'EdgeColor','r');
%% 映射到缩小的图像上
uu = s(r).BoundingBox+[0,h*0.2,0,0];%%%映射到原图
% uu = uu*2;
rectangle('Position',uu,'EdgeColor','r');
% set(gca,'units','pixels','Visible','off');
frame=getframe(gcf,[83,58,w,h]);
im=frame2im(frame);
filename = [ 'video\back' num2str(ii) '.png' ];
imwrite(im,filename);
end
if flag ==0
frame=getframe(gcf,[83,58,w,h]);
im=frame2im(frame);
filename = [ 'video\back' num2str(ii) '.png' ];
imwrite(im,filename);
end
end
评论0
最新资源