##由于作者比较菜就使用了别人的代码做了点改进和实验
原始视频(bmp序列):(猴子和狗进行了友好的握手之后,猴子对狗挥手道别,狗目送猴离去)
难点:1.画质差 2.镜头晃动 3.背景复杂
首先,使用的是https://www.cnblogs.com/tiandsp/archive/2012/07/16/2593883.html里使用的HSoptflow.m函数:
function [us,vs] = HSoptflow(Xrgb,n) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Author: Gregory Power [email protected] % 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; %Calculate gradients 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); I=imread(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
很混乱,而且尺寸很小,不会处理,于是采用了第二种方法:
原文链接:https://blog.csdn.net/garfielder007/article/details/50441756
其中缺少的computeColor.m链接:http://www.codeforge.cn/read/215288/computeColor.m__html
function img = computeColor(u,v) % computeColor color codes flow field U, V % According to the c++ source code of Daniel Scharstein % Contact: [email protected] % Author: Deqing Sun, Department of Computer Science, Brown University % Contact: [email protected] % $Date: 2007-10-31 21:20:30 (Wed, 31 Oct 2006) $ % Copyright 2007, Deqing Sun. % % All Rights Reserved % % 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); rad = sqrt(u.^2+v.^2); 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; idx = rad <= 1; col(idx) = 1-rad(idx).*(1-col(idx)); % increase saturation with radius 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 % http://members.shaw.ca/quadibloc/other/colint.htm 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;
并使用了基于局部距离的图像显著性做遮罩,即图像显著性高的地方为1,显著性低的为0。对光流法所得的显著性做遮罩,只保留图像显著性为1的地方。
最后结果如下:
实现代码:
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); Iref = imread(name1); Iinput = imread(name2); 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); BW = boundarymask(label); %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.Method = 'Lucas-Kanade'; 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=imread(name1); I2=imread(name2); 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)); v=imadjust(v,[],[],1.1);%增加对比度 亮的更亮 暗的更暗 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; img2=imadjust(img2,[],[],2); v2=imadjust(v2,[],[],2); 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
其中求显著性图像的代码为suptype.m
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); BW = boundarymask(label); %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;
并使用了该论文中的部分函数:
论文代码下载地址:http://people.csail.mit.edu/celiu/OpticalFlow/OpticalFlow.zip
最后使用ae伪造一个效果图(没错这其实是一个如何用ae伪造结果的博客):
第一步,导入bmp序列:
第二步:把序列导入合成,并点 效果->抠像->颜色差值键
第三步:选择效果控件中 视图->已校正遮罩 并使用预览图中的三个吸管(用第二个取背景 第三个取前景)求出最后效果如下:
多取几次,使猴子基本为黑,背景基本为白。
第四步:用矩形工具建立如下大小左右的遮罩。
第五步:根据猴子的移动设置对遮罩关键帧“跟踪”即可
简单做几个关键帧即可,使得猴子全身基本在框内就可以。最后直接渲染即可。
结果图:
(是不是比前边几个效果好太多 [/斜眼笑] )