本文介绍了一种基于颜色直方图的粒子滤波目标跟踪方法，详细展示了MATLAB代码实现过程，适用于视频中的运动目标估计。
视频运动目标估计是一个非线性、非高斯的过程，而粒子滤波目标跟踪由于其对于非线性和非高斯的独特特点，已被广泛应用于视频目标跟踪领域。下面主要是基于颜色直方图的粒子滤波目标跟踪结果和程序，颜色通道为hsv通道，量化bin为16 ：16:16。
下面是代码：
main.m文件代码：
function main()
mov=VideoReader('F:\两个全景\两个全景2018.4.2\剪辑好的\教十\教十350.avi'); % 更换视频目录即可
Startframe=3000;
% -----------------------------------粒子滤波变量初始化-------------------------------------
N=99;
num=1; % num为要跟踪的目标个数 , 粒子个数
delta=0.7; % 状态转移矩阵
A = [1 delta 0 0 0 0 ; 0 1 0 0 0 0 ; 0 0 1 delta 0 0 ; 0 0 0 1 0 0 ; 0 0 0 0 1 0 ; 0 0 0 0 0 1];
B = [delta^2/2 0 0 0 0 0;0 delta 0 0 0 0;0 0 delta^2/2 0 0 0;0 0 0 delta 0 0;0 0 0 0 1 0;0 0 0 0 0 1];
%
img1=read(mov,Startframe);
% 手动选取目标，并输出目标框的中心坐标以及包含目标框宽高的变量 图1变量初始化
[temp1,rect1,center1,loc1]=Get_Target1(img1,num); % loc1存放的是图1中的目标框的中心点坐标以及框的宽高
[X01,XX1,XX01,XX0_temp1,x_temp1,hist1]=state_initial1(loc1,temp1,N,num);% 粒子状态初始化%% 读取序列图像
% 初始化绘制结果的矩阵
Position1=cell(1,num);
Estmation_x1=cell(1,num);
Estmation_y1=cell(1,num);
w1=cell(1,num);
for i=1:num
% 初始化权值矩阵
w1{1,i}=ones(1,N);
end
save_hist1=cell(1,3); % 这个数组为后面处理遮挡后更新hist用
for frame=(Startframe+1):(Startframe+500)
Im1 =read(mov, frame);
for k=1:num
% 跟踪处理
[XX1{1,k},value1(k),w1{1,k},Estmation_x1{1,k},Estmation_y1{1,k},loc1{1,k}]=TrackingProcessing1(Im1,hist1{1,k},A,B,XX1{1,k},loc1{1,k},N,w1{1,k});
%状态更新
[hist1{1,k},XX01{1,k},x_temp1{1,k},Position1{1,k},save_hist1{1,k}]=UPdate1(Im1,value1(k),Estmation_x1{1,k},Estmation_y1{1,k},loc1{1,k},rect1{1,k},frame,Startframe,XX01{1,k},x_temp1{1,k},hist1{1,k},save_hist1{1,k});
%Resample
[XX1{1,k},w1{1,k}]=resample_particles1(XX1{1,k},w1{1,k},N);
% 重采样后粒子权重归一化
% w1{1,k}=ones(1,N);
end
% 显示结果
figure(1);
imshow(Im1);set(gca, 'position', [0 0 1 1 ]);axis normal;
DrawResult1(num,Position1,Estmation_x1,Estmation_y1,XX1,frame+1);
% title('Camera1');
% AFrame=getframe(gcf);
% %存储调整过大小的图片
% imwrite(AFrame.cdata,strcat('F:\粒子滤波相关\HSV求直方图的\视频的\HSV求直方图的\2017.5.18备份\1\',num2str(frame),'.jpg'))
end
function [temp,rect,center,location]=Get_Target1(rgb1,num)
%------这个函数文件是针对于读取图片的-------
% 手动选取目标，并输出目标框的中心坐标以及包含目标框
% 宽高的变量
temp=cell(1,num);
rect=cell(1,num);
center=cell(1,num);
location=cell(1,num);
for kk=1:num
figure(1),
imshow(rgb1);
[temp{1,kk},rect{1,kk}]=imcrop(rgb1);
% 得到目标框中西的坐标
center{1,kk}=[rect{1,kk}(1)+rect{1,kk}(3)/2,rect{1,kk}(2)+rect{1,kk}(4)/2];
location{1,kk}=[center{1,kk}(1),center{1,kk}(2),rect{1,kk}(3),rect{1,kk}(4)];
end
close(gcf);
end
function [X0,XX,XX0,XX0_temp,x_temp,hist1]=state_initial1(location,temp,N,num)
% 粒子状态初始化
X0=cell(1,num);
XX=cell(1,num);
XX0=cell(1,num);
XX0_temp=cell(1,num);
x_temp=cell(1,num);
hist1=cell(1,num);
for ii=1:num
X0{1,ii}=[location{1,ii}(1),0,location{1,ii}(2),0,size(temp{1,ii},1),size(temp{1,ii},2)];
XX{1,ii}=zeros(6,N);
for k=1:N
XX{1,ii}(:,k)=X0{1,ii}';
end
XX0{1,ii}=X0{1,ii};
XX0_temp{1,ii}=zeros(1,6);
x_temp{1,ii}=X0{1,ii};
%%计算权值矩阵
[hist1{1,ii},m_wei1{1,ii}]=calc_hist1(temp{1,ii});
end
end
function [XX,value,w,Estmation_x,Estmation_y,location]=TrackingProcessing1(Im1,hist1,A,B,XX,location,N,w)
% Im1 要跟踪的图像
% hist1 图像中目标的直方图
% A 状态转移矩阵
% B 转移矩阵
% XX 所有粒子初始状态
% location 目标框中心位置坐标及目标框的宽高
% N 粒子个数
% w 粒子权重
[HistDist,XX,w]=Cale_HistDist1(Im1,hist1,A,B,XX,location,N,w); % 计算每个粒子与目标框的相似度以及权重
% 得到相似度最大的粒子的位置及最大相似度值
w=w/sum(w); % 权值归一化
[Estmation_x,Estmation_y]=estimation1(XX,w); % 加权估计目标中心位置
[Estmation_x,Estmation_y]=BounderProcess1(Estmation_x,Estmation_y,location,Im1); % 防止出界
Position=[Estmation_x-location(3)/2,Estmation_y-location(4)/2,location(3),location(4)];
temp1=imcrop(Im1,round(Position));
[hist_temp1,~]=calc_hist1(temp1);
% 计算两个直方图的相似度
value=Cal_D_Between_Hist(hist_temp1,hist1);
end
function [HistDist,XX,w]=Cale_HistDist1(Im,hist1,A,B,XX,location,N,w)
% 该函数实现计算各个粒子的权值以及与上一帧样本的相似度
% 输入： Im 输入的图像
% hist1 初始的目标框的直方图向量
% A 状态转移矩阵
% B 加入的高斯扰动
% X0 上一帧的粒子的状态向量
% rect 上一帧保存的最佳粒子的位置坐标目标框的宽度
% N 粒子个数
% w 上一帧各个粒子的权重
% 输出：
% HistDist 单独用巴氏距离计算的各个粒子与上一帧目标框的相似度大小
% XX 计算得到的各个粒子的状态信息
% w 重新计算的各个粒子的权值
% 变量初始化
HistDist=zeros(1,N);
d=zeros(1,N);
p=zeros(1,N);
w_new=zeros(1,N);
% 计算图像颜色方差
% [sigma,bars]=histogram(Im);
% 噪声方差
sigma=sqrt(0.01);
for ii=1:N
W=randn(6,1);
% 将生成的N个粒子依次代入状态转移方程
XX(:,ii)=A*XX(:,ii)+15*B*W;
% 防止访问图像越界
if XX(1,ii)<1
XX(1,ii)=1;
end
if XX(1,ii)>size(Im,2)
XX(1,ii)=size(Im,2);
end
if XX(3,ii)<1
XX(3,ii)=1;
end
if XX(3,ii)>size(Im,1)
XX(3,ii)=size(Im,1);
end
location(1)=XX(1,ii);
location(2)=XX(3,ii);
location=round(location);
rect=[location(1)-location(3)/2,location(2)-location(4)/2,location(3),location(4)];
temp1=imcrop(Im,rect);
[hist2,~]=calc_hist1(temp1);
Sum1=sum(hist1);
Sum2=sum(hist2);
Sumup = sqrt(hist1.*hist2);
SumDown = sqrt(Sum1*Sum2);
Sumup = sum(Sumup);
% 由于巴士距离求出的值越小表示相似度越高，这里进行修改，用1一减，值越大则相似性越高
HistDist(ii)=1-sqrt(1-Sumup/SumDown);
% 下面的是利用课本上的求解粒子的权值
d(ii)=sqrt(1-Sumup/SumDown);
p(ii)=(1/sqrt(2*pi*sigma^2))*exp((-1*d(ii)^2)/(2*pi*sigma^2));
w_new(ii)=w(ii)*p(ii);
end
w=w_new;
end
function [ Estmation_x,Estmation_y]=estimation1(XX,w)
x_new=XX(1,:);
y_new=XX(3,:);
%利用加权平均权值估计目标新位置
Estmation_x=x_new*w';
Estmation_y=y_new*w';
end
function [Estmation_x,Estmation_y]=BounderProcess1(Estmation_x,Estmation_y,location,Im)
% 边界处理，防止粒子状态更新得到的新位置在图像外面
if Estmation_x-location(3)/2<1
Estmation_x=location(3)/2+1;
end
if Estmation_x-location(3)/2>size(Im,2)
Estmation_x=location(3)/2+size(Im,2);
end
if Estmation_y-location(4)/2<1
Estmation_y=location(4)/2+1;
end
if Estmation_y-location(4)/2>size(Im,1)
Estmation_y=location(4)/2+size(Im,1);
end
end
function value=Cal_D_Between_Hist(hist1,hist2)
Sum1=sum(hist1);
Sum2=sum(hist2);
Sumup = sqrt(hist1.*hist2);
SumDown = sqrt(Sum1*Sum2);
Sumup = sum(Sumup);
value=1-sqrt(1-Sumup/SumDown);
end
function [hist,k_m_wei] = calc_hist1(temp)
% 转为hsv通道
temp=round(255*(rgb2hsv_mex(double(temp))));
%hist为绘制直方图
[a,b,c]=size(temp);
y(1)=a/2;
y(2)=b/2;
k_m_wei=zeros(a,b);%权值矩阵
h=y(1)^2+y(2)^2 ;%带宽
%%计算权值矩阵
for i=1:a
for j=1:b
dist=(i-y(1))^2+(j-y(2))^2;
m_wei=dist/h;
if m_wei<=1
k_m_wei(i,j)=1-m_wei*m_wei; % 此种颜色直方图的计算是参考 基于分块粒子滤波的目标跟踪算法
else % 的研究 刘晓龙的硕士论文中的
k_m_wei(i,j)=0;
end
end
end
C=1/sum(sum(k_m_wei)); % 归一化系数
%%计算目标权值直方图
hist=zeros(1,4096);
for i=1:a
for j=1:b %rgb颜色空间量化为16*16*16bins
q_h=floor(double(temp(i,j,1))/16);%fix是趋于0取整函数
q_s=floor(double(temp(i,j,2))/16);
q_v=floor(double(temp(i,j,3))/16);
q_temp=q_h*16+q_s*256+1*q_v;%设置每个像素点红色、绿色、蓝色分量所占的比重
hist(q_temp+1)= hist(q_temp+1)+k_m_wei(i,j);%计算直方图中每个像素点所占的权重
%(q_temp+1)=hist(q_temp+1)*C;
end
end
hist=hist*C;
end
function [hist1,XX0,x_temp,Position1,save_hist]=UPdate1(Im,value,Estmation_x,Estmation_y,location,rect,frame,Startframe,XX0,x_temp,hist,save_hist)
% 该函数实现对粒子状态的更新
% 输入： value 主程序里面求出的相似度最大的值
% Estmation_x 表示估计当前帧的粒子的目标框中心位置横坐标
% Estmation_y 表示估计当前帧的粒子的目标框中心位置纵坐标
% XX 主程序里面求出的每个粒子的状态矩阵
% XX0 主程序里面初始的最初粒子状态
% XX0_temp 6*1的全零矩阵
% x_temp 主程序里面初始的最初粒子状态
% frame 当前帧的图像
% startframe 当前帧所在帧数
% location 存放抠取的图像的宽高即坐标的数组
% rect 主函数里面存放抠取的图像的宽高即坐标的数组
% 输出： vv 这个结构体保存的是画矩形框所用的左上角的坐标及矩形框的宽高
% X0 初始状态的更新
% XX0 初始状态的更新
XX0_temp=zeros(1,6);
if value>0.8
% 如果相似度大于0.7，直接用估计的目标位置作为下一帧的初始目标
rect(1)=Estmation_x-location(3)/2;
rect(2)=Estmation_y-location(4)/2;
temp1=imcrop(Im,round(rect));
[hist1,~]=calc_hist1(temp1);
XX0(1)=Estmation_x;
XX0(3)=Estmation_y;
XX0_temp(2)=( XX0_temp(2)+XX0(1)-x_temp(1) );
XX0_temp(4)=( XX0_temp(4)+XX0(3)-x_temp(3) );
XX0(2)=XX0_temp(2)/(frame-Startframe);
XX0(4)=XX0_temp(4)/(frame-Startframe);
x_temp=XX0;
Position1=[rect(1),rect(2),location(3),location(4)];
else
% 否则结合上一帧目标的位置以及速度估计下一帧的初始位置
XX0(1)= x_temp(1)+x_temp(2); % 上一帧粒子的位置加上粒子的速度（即一帧的位移长度）
XX0(3)= x_temp(3)+x_temp(4);
XX0_temp(2)=(XX0_temp(2)+XX0(1)-x_temp(1)); % 重新估计粒子的速度
XX0_temp(4)=(XX0_temp(4)+XX0(3)-x_temp(3));
XX0(2)=XX0_temp(2)/(frame-Startframe);
XX0(4)=XX0_temp(4)/(frame-Startframe);
x_temp=XX0;
% 若value小于阈值 判断目标已经被遮挡 输出下面这个
% 输出目标位置
Position1(1,1)=Estmation_x-location(3)/2;
Position1(1,2)=Estmation_y-location(4)/2;
Position1(1,3)=location(3);
Position1(1,4)=location(4);
% 得到目标位置对应的区域的直方图
temp1=imcrop(Im,round(Position1));
[hist_temp1,~]=calc_hist1(temp1);
save_hist=[save_hist;hist_temp1];
% 对hist1直接赋值遮挡前的直方图
hist1=hist;
% 如果满足下面条件 更改hist1赋值
% 本程序里面更新hist的条件是：
% 如果判断了目标已经被遮挡，保存被遮挡前的直方图，然后依次保存连续三帧的
% 估计的目标区域的直方图，若连续两帧的直方图的相似度大于0.7并且第三个估计
% 的直方图与遮挡时保存的直方图的相似度大于0.5，即判断目标已经跟踪回来，用
% 第三个估计的直方图代替之前保存的遮挡时保存的直方图
if size((save_hist),1)==3
tmp1=Cal_D_Between_Hist(save_hist(1,:),save_hist(2,:));
tmp2=Cal_D_Between_Hist(save_hist(2,:),save_hist(3,:));
tmp3=Cal_D_Between_Hist(save_hist(3,:),hist);
if tmp1>0.75&tmp2>0.75&tmp3>0.6
hist1=save_hist(3,:);
end
save_hist=[];
end
end
end
function [XX_resample,w_new]=resample_particles1(XX,w,N)
% 该函数实现对粒子的重采样
% 输入： w 粒子权重
% XX 主程序里面求出的每个粒子的状态矩阵
% N 粒子个数
% 输出： w 粒子权重
% XX 主程序里面求出的每个粒子的状态矩阵
XX_resample=XX;
w_new=w;
% 判断是否需要重采样
% Neff表示有效的粒子个数(也即是否需要重采样的阈值)
Neff=1/sum(w.*w);
if Neff<50
for i = 1 : N
u = rand; % uniform random number between 0 and 1
qtempsum = 0;
for j = 1 : N
qtempsum = qtempsum + w(j);
if qtempsum >= u
XX_resample(:,i)=XX(:,j);
break;
end
end
end
w_new=ones(1,N)/N;
end
end
function DrawResult1(num,Position1,Estmation_x,Estmation_y,XX,frame)
% function DrawResult1(num,Position1,Estmation_x,Estmation_y,XX,frame,EpipolarLine)
% 定义不同颜色
Color=[1 0 0;
0 1 0;
0 0 1
1 1 0;
1 0 1;
0 1 1;
1 1 1];
for k=1:num
hold on;
rectangle('position',Position1{1,k},'EdgeColor',Color(k,:));
hold on;
plot(XX{1,k}(1,:),XX{1,k}(3,:),'r.')
hold on;
plot(Estmation_x{1,k},Estmation_y{1,k},'g+','MarkerSize',10);
end
text(10,20,int2str(frame),'FontSize',15,'Color','r');
end
endrgb2hsv_mex.c文件代码。
/*
Convert an image from RGB to HSV
Usage: A = rgb2hsv_mex(im)
-----
Inputs
im Image (n x m x 3) in double format
Ouputs
Z HSV image (n x m x 3)
Example
im = double(ceil(255*rand(200 , 300 , 3)));
Z = rgb2hsv_mex(im);
To compile
mex rgb2hsv_mex.c
mex -f mexopts_intel10.bat -output rgb2hsv_mex.dll rgb2hsv_mex.c
-------
*/
#include <math.h>
#include "mex.h"
#ifndef MAX
#define MIN(a , b) ((a) < (b) ? (a) : (b))
#define MAX(a , b) ((a) >= (b) ? (a) : (b))
#endif
#define NO_HUE 0
/*------------------------ Headers --------------------------------------*/
void rgb2hsv (double * , double * , int , int);
/*-------------------------------------------------------------------------*/
void mexFunction( int nlhs, mxArray *plhs[] , int nrhs, const mxArray *prhs[] )
{
double *RGB;
double *HSV;
const int *dimsRGB;
int n , m , numdimsRGB;
/*---------------------------------------------------------------*/
/*---------------------- PARSE INPUT ----------------------------*/
/*---------------------------------------------------------------*/
if(nrhs != 1)
{
mexErrMsgTxt("1 input is required");
}
/*---------------- Input 1 ----------------------*/
RGB = mxGetPr(prhs[0]);
numdimsRGB = mxGetNumberOfDimensions(prhs[0]);
dimsRGB = mxGetDimensions(prhs[0]);
if (numdimsRGB != 3)
{
mexErrMsgTxt("First input must at least a Matrix (n x m x 3)");
}
n = dimsRGB[0];
m = dimsRGB[1];
/*------------------- Output 1 ----------------------*/
plhs[0] = mxCreateNumericArray(numdimsRGB, dimsRGB, mxDOUBLE_CLASS, mxREAL);
HSV = mxGetPr(plhs[0]);
/*----------- Convert into HSV ---------------*/
rgb2hsv(RGB , HSV , n , m);
/*--------------------------------------------*/
}
/* ----------------------------------------------------------------------------------------------- */
void rgb2hsv (double *RGB , double *HSV , int n , int m)
{
double r , g , b , max , min , delta;
double h , s , v , cte = 1.0/255.0, cte1 = 1.0/360.0;
int nm = n*m , i;
for (i = 0 ; i < nm ; i++)
{
r = RGB[i]*cte;
g = RGB[i + nm]*cte;
b = RGB[i + 2*nm]*cte;
max = MAX(r , MAX(g , b));
min = MIN(r , MIN(g , b));
delta = (max - min);
v = max;
if (max != 0.0)
{
s = delta/max;
}
else
{
s = 0.0;
}
if (s != 0.0)
{
if (r == max)
{
h = (g - b)/delta;
}
else if (g == max)
{
h = 2.0 + (b - r)/delta;
}
else if (b == max)
{
h = 4.0 + (r - g)/delta;
}
h *= 60.0;
if (h < 0)
{
h += 360.0;
}
h *= cte1;
}
else
{
h = NO_HUE;
}
HSV[i] = h;
HSV[i + nm] = s;
HSV[i + 2*nm] = v;
}
}
/* ----------------------------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------------------------------- */
以上是所有代码。包括两个部分，一个是主函数main.m，一个是计算rgb通道转hsv通道的c代码。运行上面代码之前，先将将上面代码保存的m文件和该.c文件放置在同一个文件夹之中，然后请确保你的MATLAB可以编译.c文件，具体如何编译，请找度娘。用MATLAB对.c文件进行编译，编译完之后，在你的当前文件夹会生成一个rgb2hsv_mex.mexw64或rgb2hsv_mex.mexw32的文件。具体是哪个，由你的电脑位数决定。然后修改main.m之中的视频文件路径，运行即可。（该代码是针对视频文件进行处理的，https://download.youkuaiyun.com/download/u011624019/10140469，该网址是所有打包下载，但是针对于64位电脑才可以运行的）
扫一扫
362
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
