利用ICA 进行运动跟踪

利用ICA的原理，可以分离两幅图像中不相同的部分，这就为运动跟踪，建立了条件。

下面通过一个案例进行说明，代码采用matlab。

Matlab代码如下：

videoDemo.m

video = mmreader('1.avi');%% 读取要处理的视频nFrames = video.NumberOfFrames;   %得到帧数H = video.Height;     %得到高度W = video.Width;      %得到宽度Rate = video.FrameRate; %得到视频帧率% Preallocate movie structure.mov(1:nFrames) = struct('cdata',zeros(H,W,3,'uint8'),'colormap',[]);back = rgb2gray(imread('back.png'));%% 读取背景back = imresize(back,0.5);%% 缩小原始图像%read one frame every timefor i = 222:2:3100   mov(i).cdata = read(video,i);    P= mov(i).cdata;    P =rgb2gray(P);%% 彩色图像转换到灰度图像    d = imresize(P,0.5);%     disp('当前播帧数：'),disp(i);    dst = compute(d,back,P,i);%%  视频处理            pause(0.05);%%  停顿视频end

其中，compute.m

function [ dst ] = compute( I1,I2 ,src,ii )%UNTITLED3 Summary of this function goes here%   Detailed explanation goes heremaxV1 = double(max(max(I1)))/255;%%获得原始图像的最大最小值minV1 = double(min(min(I1)))/255;%%<span style="font-family: Arial, Helvetica, sans-serif;">获得原始图像的最大最小值</span>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);    endend%%%%%%%%%%%%%%%%%%%%%%%%%%%  白化  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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;                                % 保存所得向量bend%%%%%%%%%%%%%%%%%%%%%%%%%%  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<minV1dst1 = imadjust(reshape(ICAedS(1,:),len1,len2),[minv1,maxv1],[minV1,maxV1],1);elsedst1 = imadjust(reshape(ICAedS(1,:),len1,len2));endif minv2<minV2dst2 = imadjust(reshape(ICAedS(2,:),len1,len2),[minv2,maxv2],[minV2,maxV2],1);elsedst2 = 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);               endif flag ==0    frame=getframe(gcf,[83,58,w,h]);     im=frame2im(frame);     filename = [  'video\back' num2str(ii) '.png' ];        imwrite(im,filename);endend

视频的处理效果如下图所示：只判断车所在的那条路上的情况，简单的demo程序在最后面的链接中。

我将工程的所有内容放到http://download.csdn.net/detail/yeyang911/9227303 这个链接上，感兴趣的可以下载运行。