# 使用光流法进行视频显著性检测（Matlab）(上)

##由于作者比较菜就使用了别人的代码做了点改进和实验

function [us,vs] = HSoptflow(Xrgb,n)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Gregory Power gregory.power@wpafb.af.mil
% This MATLAB code shows a Motion Estimation map created by
% using a Horn and Schunck motion estimation technique on two
% consecutive frames.  Input requires.
%     Xrgb(h,w,d,N) where X is a frame sequence of a certain
%                height(h), width (w), depth (d=3 for color),
%                and number of frames (N).
%     n= is the starting frame number which is less than N
%     V= the output variable which is a 2D velocity array
%
% Sample Call: V=HSoptflow(X,3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[h,w,d,N]=size(Xrgb);
if n>N-1
error(1,'requested file greater than frame number required');
end;
%get two image frames
if d==1
Xn=double(Xrgb(:,:,1,n));
Xnp1=double(Xrgb(:,:,1,n+1));
elseif d==3
Xn=double(Xrgb(:,:,1,n)*0.299+Xrgb(:,:,2,n)*0.587+Xrgb(:,:,3,n)*0.114);
Xnp1=double(Xrgb(:,:,1,n+1)*0.299+Xrgb(:,:,2,n+1)*0.587+Xrgb(:,:,3,n+1)*0.114);
else
error(2,'not an RGB or Monochrome image file');
end;

%get image size and adjust for border
size(Xn);
hm5=h-5; wm5=w-5;
z=zeros(h,w); v1=z; v2=z;

%initialize
dsx2=v1; dsx1=v1; dst=v1;
alpha2=625;
imax=20;

dst(5:hm5,5:wm5) = ( Xnp1(6:hm5+1,6:wm5+1)-Xn(6:hm5+1,6:wm5+1) + Xnp1(6:hm5+1,5:wm5)-Xn(6:hm5+1,5:wm5)+ Xnp1(5:hm5,6:wm5+1)-Xn(5:hm5,6:wm5+1) +Xnp1(5:hm5,5:wm5)-Xn(5:hm5,5:wm5))/4;
dsx2(5:hm5,5:wm5) = ( Xnp1(6:hm5+1,6:wm5+1)-Xnp1(5:hm5,6:wm5+1) + Xnp1(6:hm5+1,5:wm5)-Xnp1(5:hm5,5:wm5)+ Xn(6:hm5+1,6:wm5+1)-Xn(5:hm5,6:wm5+1) +Xn(6:hm5+1,5:wm5)-Xn(5:hm5,5:wm5))/4;
dsx1(5:hm5,5:wm5) = ( Xnp1(6:hm5+1,6:wm5+1)-Xnp1(6:hm5+1,5:wm5) + Xnp1(5:hm5,6:wm5+1)-Xnp1(5:hm5,5:wm5)+ Xn(6:hm5+1,6:wm5+1)-Xn(6:hm5+1,5:wm5) +Xn(5:hm5,6:wm5+1)-Xn(5:hm5,5:wm5))/4;

for i=1:imax
delta=(dsx1.*v1+dsx2.*v2+dst)./(alpha2+dsx1.^2+dsx2.^2);
v1=v1-dsx1.*delta;
v2=v2-dsx2.*delta;
end
u=z; u(5:hm5,5:wm5)=v1(5:hm5,5:wm5);
v=z; v(5:hm5,5:wm5)=v2(5:hm5,5:wm5);

xskip=round(h/64);
[hs,ws]=size(u(1:xskip:h,1:xskip:w));
us=zeros(hs,ws); vs=us;

N=xskip^2;
for i=1:hs-1
for j=1:ws-1
hk=i*xskip-xskip+1;
hl=i*xskip;
wk=j*xskip-xskip+1;
wl=j*xskip;
us(i,j)=sum(sum(u(hk:hl,wk:wl)))/N;
vs(i,j)=sum(sum(v(hk:hl,wk:wl)))/N;
end
end

figure(1);
quiver(us,vs);
colormap('default');
axis ij;
axis tight;
axis equal;

close all;
clear all;
clc;
img_dir = dir('./monkeydog/*.bmp');
[m n]=size(img_dir);
raw=zeros(240,320,1,m);
for i=1:m
name=strcat('./monkeydog/',img_dir(i).name);
raw(:,:,1,i)=rgb2gray(I);
end

for i=1:m-1
V=HSoptflow(raw,i);
imshow(V)
imgname=strcat('C:\Users\71405\Desktop\回收站\计算机视觉\2018.6.12\output2\img_',num2str(i),'.jpg');
saveas(gcf,imgname);#保存图片
end

function img = computeColor(u,v)

%   computeColor color codes flow field U, V

%   According to the c++ source code of Daniel Scharstein
%   Contact: schar@middlebury.edu

%   Author: Deqing Sun, Department of Computer Science, Brown University
%   Contact: dqsun@cs.brown.edu
%   $Date: 2007-10-31 21:20:30 (Wed, 31 Oct 2006)$

%
%
% Permission to use, copy, modify, and distribute this software and its
% documentation for any purpose other than its incorporation into a
% commercial product is hereby granted without fee, provided that the
% above copyright notice appear in all copies and that both that
% copyright notice and this permission notice appear in supporting
% documentation, and that the name of the author and Brown University not be used in
% advertising or publicity pertaining to distribution of the software
% without specific, written prior permission.
%
% THE AUTHOR AND BROWN UNIVERSITY DISCLAIM ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
% INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ANY
% PARTICULAR PURPOSE.  IN NO EVENT SHALL THE AUTHOR OR BROWN UNIVERSITY BE LIABLE FOR
% ANY SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
% WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN AN
% ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF
% OR IN CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE.

nanIdx = isnan(u) | isnan(v);
u(nanIdx) = 0;
v(nanIdx) = 0;

colorwheel = makeColorwheel();
ncols = size(colorwheel, 1);

a = atan2(-v, -u)/pi;

fk = (a+1) /2 * (ncols-1) + 1;  % -1~1 maped to 1~ncols

k0 = floor(fk);                 % 1, 2, ..., ncols

k1 = k0+1;
k1(k1==ncols+1) = 1;

f = fk - k0;

for i = 1:size(colorwheel,2)
tmp = colorwheel(:,i);
col0 = tmp(k0)/255;
col1 = tmp(k1)/255;
col = (1-f).*col0 + f.*col1;

col(~idx) = col(~idx)*0.75;             % out of range

img(:,:, i) = uint8(floor(255*col.*(1-nanIdx)));
end;

%%
function colorwheel = makeColorwheel()

%   color encoding scheme

%   adapted from the color circle idea described at

RY = 15;
YG = 6;
GC = 4;
CB = 11;
BM = 13;
MR = 6;

ncols = RY + YG + GC + CB + BM + MR;

colorwheel = zeros(ncols, 3); % r g b

col = 0;
%RY
colorwheel(1:RY, 1) = 255;
colorwheel(1:RY, 2) = floor(255*(0:RY-1)/RY)';
col = col+RY;

%YG
colorwheel(col+(1:YG), 1) = 255 - floor(255*(0:YG-1)/YG)';
colorwheel(col+(1:YG), 2) = 255;
col = col+YG;

%GC
colorwheel(col+(1:GC), 2) = 255;
colorwheel(col+(1:GC), 3) = floor(255*(0:GC-1)/GC)';
col = col+GC;

%CB
colorwheel(col+(1:CB), 2) = 255 - floor(255*(0:CB-1)/CB)';
colorwheel(col+(1:CB), 3) = 255;
col = col+CB;

%BM
colorwheel(col+(1:BM), 3) = 255;
colorwheel(col+(1:BM), 1) = floor(255*(0:BM-1)/BM)';
col = col+BM;

%MR
colorwheel(col+(1:MR), 3) = 255 - floor(255*(0:MR-1)/MR)';
colorwheel(col+(1:MR), 1) = 255;

img_dir = dir('./monkeydog/*.bmp');

[m n]=size(img_dir);

for q=1:m
name1=strcat('./monkeydog/',img_dir(q).name);
name2=strcat('./monkeydog/',img_dir(q+1).name);
I=Iinput;

%% 求后一帧的显著性
[m n d]=size(I);
Iq=imresize(uint8(I),[240,240*n/m]);
Iref=imresize(uint8(Iref),[240,240*n/m]);
Iinput=imresize(uint8(Iinput),[240,240*n/m]);
mode=fspecial('gaussian', 1, 1);
%IS=imfilter(Iq,mode,'replicate');
IS=Iq;
[label,num]=superpixels(IS,300);
%imshow(imoverlay(IS,BW,'cyan'),'InitialMagnification',100);
%figure
[m n]=size(label);
supmean=uint64(zeros(num,5));
labelnum=zeros(num);

for i=1:m
for j=1:n
supmean(label(i,j),1)=supmean(label(i,j),1)+i;
supmean(label(i,j),2)=supmean(label(i,j),2)+j;
supmean(label(i,j),3)=(supmean(label(i,j),3)+uint64(IS(i,j,1)));
supmean(label(i,j),4)=(supmean(label(i,j),4)+uint64(IS(i,j,2)));
supmean(label(i,j),5)=(supmean(label(i,j),5)+uint64(IS(i,j,3)));
labelnum(label(i,j))=labelnum(label(i,j))+1;
end
end
for i=1:num
supmean(i,1)=uint16(supmean(i,1)./labelnum(i));
supmean(i,2)=uint16(supmean(i,2)./labelnum(i));
supmean(i,3)=uint16(supmean(i,3)./labelnum(i));
supmean(i,4)=uint16(supmean(i,4)./labelnum(i));
supmean(i,5)=uint16(supmean(i,5)./labelnum(i));
end

IM=zeros(m,n,3);
for i=1:m
for j=1:n
IM(i,j,1)=uint8(supmean(label(i,j),3));
IM(i,j,2)=uint8(supmean(label(i,j),4));
IM(i,j,3)=uint8(supmean(label(i,j),5));
end
end
supim=uint8(cat(3,IM(:,:,1),IM(:,:,2),IM(:,:,3)));
%imshow(supim);

%全局
supmean=double(supmean);
dist_all=zeros(num,1);
w0=0.0;
for i=1:num
for j=1:num
if ((supmean(i,1)-supmean(j,1))^2+(supmean(i,2)-supmean(j,2))^2)<500
dist_all(i)=dist_all(i)+(w0*(supmean(i,1)-supmean(j,1))^2+w0*(supmean(i,2)-supmean(j,2))^2+(supmean(i,3)-supmean(j,3))^2+(supmean(i,4)-supmean(j,4))^2+(supmean(i,5)-supmean(j,5))^2)^0.5;
end
end
end
%归一化
dist_min=min(dist_all);
dist_max=max(dist_all);
for i=1:num
dist_all(i)=(dist_all(i)-dist_min)*255/(dist_max-dist_min);
end
dist_all=uint8(dist_all);
im_all=zeros(240,240*n/m);
for i=1:m
for j=1:n
im_all(i,j)=dist_all(label(i,j));
end
end
im_all=uint8(im_all);
im_all_bw=im2bw(im_all,0.55);
%im_all=im2bw(im_all,0.05);

%% 光流法
if size(Iref, 3) == 3
Irefg = rgb2gray(Iref);
Iinputg = rgb2gray(Iinput);
else
Irefg = Iref;
Iinputg = Iinput;
end
% 创建光流对象及类型转化对象
opticalFlow = vision.OpticalFlow('ReferenceFrameDelay',1);
converter = vision.ImageDataTypeConverter;
% 修改光流对象的配置
opticalFlow.OutputValue = 'Horizontal and vertical components in complex form'; % 返回复数形式光流图
opticalFlow.ReferenceFrameSource = 'Input port'; % 对比两张图片，而不是视频流
if 1 % 使用的算法
opticalFlow.NoiseReductionThreshold = 0.001; % 默认是0.0039
else
opticalFlow.Method = 'Horn-Schunck';
opticalFlow.Smoothness = 0.5; % 默认是1
end
% 调用光流对象计算两张图片的光流
Iinputg_c = step(converter, Iinputg);
Irefg_c = step(converter, Irefg);
of = step(opticalFlow, Iinputg_c, Irefg_c);
ofH = real(of);
ofV = imag(of);

[m n]=size(ofV);
for j=1:m
for k=1:n
if abs(double(ofV(m,n)))<0.5
ofV(m,n)=1;
else
ofV(m,n)=0;
end
end
end
%im_all=1-im_all;
% 将光流转化为彩色图显示
% B=strel('disk',1);
%ofV=imdilate(ofV,B);
%ofV=imerode(ofV,B);
ofI = computeColor(ofH, ofV);
ofI(:,:,1)=ofI(:,:,1).*uint8(im_all_bw);
ofI(:,:,2)=ofI(:,:,2).*uint8(im_all_bw);
ofI(:,:,3)=ofI(:,:,3).*uint8(im_all_bw);
ofI=cat(3,ofI(:,:,1),ofI(:,:,2),ofI(:,:,3));
%显示结果
%     ofV=uint8(ofV);
%     ofV=ofV.*im_all;
for qa=1:m
for qb=1:n
if ofI(qa,qb,1)+ofI(qa,qb,2)+ofI(qa,qb,3)==0
ofI(qa,qb,1)=255;
ofI(qa,qb,2)=255;
ofI(qa,qb,3)=255;
end
end
end
subplot(121),imshow(I);
subplot(122),
imshow(ofI)
pause(0.01);
imgname=strcat('C:\Users\71405\Desktop\回收站\计算机视觉\2018.6.12\output3\img_',num2str(q),'.jpg');
saveas(gcf,imgname);
end


效果也不太好，于是使用Liu Ce所撰写的代码：

原文链接：http://people.csail.mit.edu/celiu/OpticalFlow/

效果图：

（左：光流为距离显著做遮罩 右：距离显著由光流做遮罩）

效果还凑活，就是有的时候会丢失，导致结果一闪一闪的。

代码：

close all;
clear all;
clc;
img_dir = dir('./monkeydog/*.bmp');
[m n]=size(img_dir);
for i=1:m-1
name1=strcat('./monkeydog/',img_dir(i).name);
name2=strcat('./monkeydog/',img_dir(i+1).name);
I1=imresize(I1,[120,160]);
I2=imresize(I2,[120,160]);
alpha = 0.012;
ratio = 0.75;
minWidth = 20;
nOuterFPIterations = 7;
nInnerFPIterations = 1;
nSORIterations = 30;
para = [alpha,ratio,minWidth,nOuterFPIterations,nInnerFPIterations,nSORIterations];
[vx,vy,warpI2] = Coarse2FineTwoFrames(I1,I2,para);
vx=vx-mean(mean(vx));%去中心化
vy=vy-mean(mean(vy));
v=sqrt(vx.^2+vy.^2);
v1=v;
[l,d]=size(v);

%v=log(double(v));
maxi=max(max(v));
mini=min(min(v));
for k=1:l
for o=1:d
v(k,o)=(v(k,o)-mini)*255/(maxi-mini);%归一化到[0 255]
end
end
v=uint8(v);
img=suptype(I2);
%方法一： v作为遮罩
v1=im2bw(v,0.85);
img3=filter2(fspecial('gaussian'),sample)
img2=0.8*img.*uint8(v1)+0.2*img;
%方法二：img作为遮罩
img1=im2bw(img,0.65);
v2=0.8*v.*uint8(img1)+0.2*img;

subplot(121),imshow(img2);
subplot(122),imshow(v2);
pause(0.01)
imgname=strcat('C:\Users\71405\Desktop\回收站\计算机视觉\2018.6.12\output\img_',num2str(i),'.jpg');
saveas(gcf,imgname);
end

function img = suptype(I)
[m n d]=size(I);
%I=imresize(uint8(I),[120,120*n/m]);
%mode=fspecial('gaussian', 1, 1);
%IS=imfilter(I,mode,'replicate');
IS=I;
[label,num]=superpixels(IS,300);
%imshow(imoverlay(IS,BW,'cyan'),'InitialMagnification',100);
%figure
[m n]=size(label);
supmean=uint64(zeros(num,5));
labelnum=zeros(num);

for i=1:m
for j=1:n
supmean(label(i,j),1)=supmean(label(i,j),1)+i;
supmean(label(i,j),2)=supmean(label(i,j),2)+j;
supmean(label(i,j),3)=(supmean(label(i,j),3)+uint64(IS(i,j,1)));
supmean(label(i,j),4)=(supmean(label(i,j),4)+uint64(IS(i,j,2)));
supmean(label(i,j),5)=(supmean(label(i,j),5)+uint64(IS(i,j,3)));
labelnum(label(i,j))=labelnum(label(i,j))+1;
end
end
for i=1:num
supmean(i,1)=uint16(supmean(i,1)./labelnum(i));
supmean(i,2)=uint16(supmean(i,2)./labelnum(i));
supmean(i,3)=uint16(supmean(i,3)./labelnum(i));
supmean(i,4)=uint16(supmean(i,4)./labelnum(i));
supmean(i,5)=uint16(supmean(i,5)./labelnum(i));
end

IM=zeros(m,n,3);
for i=1:m
for j=1:n
IM(i,j,1)=uint8(supmean(label(i,j),3));
IM(i,j,2)=uint8(supmean(label(i,j),4));
IM(i,j,3)=uint8(supmean(label(i,j),5));
end
end
supim=uint8(cat(3,IM(:,:,1),IM(:,:,2),IM(:,:,3)));
%imshow(supim);

%全局
supmean=double(supmean);
dist_all=zeros(num,1);
w0=0.0;
for i=1:num
for j=1:num
if ((supmean(i,1)-supmean(j,1))^2+(supmean(i,2)-supmean(j,2))^2)<500
dist_all(i)=dist_all(i)+(w0*(supmean(i,1)-supmean(j,1))^2+w0*(supmean(i,2)-supmean(j,2))^2+(supmean(i,3)-supmean(j,3))^2+(supmean(i,4)-supmean(j,4))^2+(supmean(i,5)-supmean(j,5))^2)^0.5;
end
end
end
%归一化
dist_min=min(dist_all);
dist_max=max(dist_all);
for i=1:num
dist_all(i)=(dist_all(i)-dist_min)*255/(dist_max-dist_min);
end
dist_all=uint8(dist_all);
im_all=zeros(120,120*n/m);
for i=1:m
for j=1:n
im_all(i,j)=dist_all(label(i,j));
end
end
im_all=uint8(im_all);
%im_all=im2bw(im_all,0.55);
%im_all=uint8(im2bw(im_all,0.55)).*I(:,:,1);

%显示图片
img=im_all;

（是不是比前边几个效果好太多 [/斜眼笑]  ）

• 广告
• 抄袭
• 版权
• 政治
• 色情
• 无意义
• 其他

120