闲着也是闲着，不如看看代码吧，最近开始接触压缩感知跟踪，是从朋友那里看到的，所以想找了一些代码跑跑看看，算是体验一下吧！！
论文：Real-Time Compressive Tracking
作者：Kaihua Zhang
单位：香港理工大学
作者主页：http://www4.comp.polyu.edu.hk/~cskhzhang/
CT主页：http://www4.comp.polyu.edu.hk/~cslzhang/CT/CT.htm
作者提供了这篇论文matlab和C++的代码（见CT主页）
运行平台： matlab2010a
下面我直接贴代码了。
mexCompile.m文件，这里利用matlab和C++的混合，这个文件将integral.cpp和FtrVal.cpp文件进行编译为.mex文件
mex integral.cpp;
mex FtrVal.cpp;Runtracker.m这是主文件，首先对初始帧进行处理，之后运行时从第二针开始进行跟踪
% Demo for paper "Real-Time Compressive Tracking,"Kaihua Zhang, Lei Zhang, and Ming-Hsuan Yang
% To appear in European Conference on Computer Vision (ECCV 2012), Florence, Italy, October, 2012
% Implemented by Kaihua Zhang, Dept.of Computing, HK PolyU.
% Email: zhkhua@gmail.com
% Date: 11/12/2011
% Revised by Kaihua Zhang, 15/12/2011
% Revised by Kaihua Zhang, 11/7/2012
clc;clear all;close all;
%----------------------------------
rand('state',0);
%----------------------------------
%添加路径
addpath('./data');
%----------------------------------
%加载文件
load init.txt;
%[x y width height]
initstate = init;%initial tracker
%----------------------------Set path
img_dir = dir('./data/*.png');
%---------------------------
%读取初始帧
img = imread(img_dir(1).name);
%将uint8转化为double类型
img = double(img(:,:,1));
%----------------------------------------------------------------
%负样本训练数目
trparams.init_negnumtrain = 50;%number of trained negative samples
%正样本搜索半径，设置为4~8
trparams.init_postrainrad = 4.0;%radical scope of positive samples
%目标初始状态[x,y,w,h]
trparams.initstate = initstate;% object position [x y width height]
%新的一帧中的搜索窗口的半径，通常设置为15~35
trparams.srchwinsz = 20;% size of search window
% Sometimes, it affects the results.
%-------------------------
% classifier parameters
%宽
clfparams.width = trparams.initstate(3);
%高
clfparams.height= trparams.initstate(4);
%-------------------------
% feature parameters
% number of rectangle from 2 to 4.
ftrparams.minNumRect = 2;
%随机矩阵的每一行中的非零入口最大值，设置为4-6
ftrparams.maxNumRect = 4;
%-------------------------
%弱分类器个数
M = 50;% number of all weaker classifiers, i.e,feature pool
%-------------------------
%均值
posx.mu = zeros(M,1);% mean of positive features
negx.mu = zeros(M,1);
%方差
posx.sig= ones(M,1);% variance of positive features
negx.sig= ones(M,1);
%-------------------------Learning rate parameter
%学习速率，通常设置为0.7~0.95
lRate = 0.85;
%-------------------------
%compute feature template
%@clfparams: 分类器参数（宽和高）
%@ftrparams：特征参数
%@M: 弱分类器个数=50
[ftr.px,ftr.py,ftr.pw,ftr.ph,ftr.pwt] = HaarFtr(clfparams,ftrparams,M);
%-------------------------
%compute sample templates
%@para: img: 图像
%@para: initstate：初始状态[x,y,w,h]
%@para: trparams.init_postrainrad正样本搜索半径，设置为4~8
%10000: 样本的最大数量
posx.sampleImage = sampleImg(img,initstate,trparams.init_postrainrad,0,100000);
negx.sampleImage = sampleImg(img,initstate,1.5*trparams.srchwinsz,4+trparams.init_postrainrad,50);
%-----------------------------------
%--------Feature extraction
%计算积分图，函数integral.mexw32
iH = integral(img);%Compute integral image
posx.feature = getFtrVal(iH,posx.sampleImage,ftr);
negx.feature = getFtrVal(iH,negx.sampleImage,ftr);
%--------------------------------------------------
[posx.mu,posx.sig,negx.mu,negx.sig] = classiferUpdate(posx,negx,posx.mu,posx.sig,negx.mu,negx.sig,lRate);% update distribution parameters
%-------------------------------------------------
num = length(img_dir);% number of frames
%--------------------------------------------------------
x = initstate(1);% x axis at the Top left corner
y = initstate(2);
w = initstate(3);% width of the rectangle
h = initstate(4);% height of the rectangle
%--------------------------------------------------------
%--------------运行时------------------------------------
%--------------------------------------------------------
for i = 2:num
%读取视频帧
img = imread(img_dir(i).name);
imgSr = img;% imgSr is used for showing tracking results.
img = double(img(:,:,1));
%计算新帧的特征模板
detectx.sampleImage = sampleImg(img,initstate,trparams.srchwinsz,0,100000);
%计算新帧的积分图
iH = integral(img);%Compute integral image
%对新帧进行计算其特征值
detectx.feature = getFtrVal(iH,detectx.sampleImage,ftr);
%------------------------------------
r = ratioClassifier(posx,negx,detectx);% compute the classifier for all samples
clf = sum(r);% linearly combine the ratio classifiers in r to the final classifier
%-------------------------------------
[c,index] = max(clf);
%--------------------------------
%更新模板状态[x,y,w,h]
x = detectx.sampleImage.sx(index);
y = detectx.sampleImage.sy(index);
w = detectx.sampleImage.sw(index);
h = detectx.sampleImage.sh(index);
initstate = [x y w h];
%-------------------------------Show the tracking results
%-----显示跟踪结果----
imshow(uint8(imgSr));
rectangle('Position',initstate,'LineWidth',4,'EdgeColor','r');
hold on;
%在视频帧中添加文字
text(5, 18, strcat('#',num2str(i)), 'Color','y', 'FontWeight','bold', 'FontSize',20);
%设置位置
set(gca,'position',[0 0 1 1]);
pause(0.00001);
hold off;
%------------------------------Extract samples
%提取样本，为下一帧做准备！！
posx.sampleImage = sampleImg(img,initstate,trparams.init_postrainrad,0,100000);
negx.sampleImage = sampleImg(img,initstate,1.5*trparams.srchwinsz,4+trparams.init_postrainrad,trparams.init_negnumtrain);
%--------------------------------------------------Update all the features
%更新特征
posx.feature = getFtrVal(iH,posx.sampleImage,ftr);
negx.feature = getFtrVal(iH,negx.sampleImage,ftr);
%--------------------------------------------------
[posx.mu,posx.sig,negx.mu,negx.sig] = classiferUpdate(posx,negx,posx.mu,posx.sig,negx.mu,negx.sig,lRate);% update distribution parameters
endHaarFtr.m文件
function [px,py,pw,ph,pwt] = HaarFtr(clfparams,ftrparams,M)
% $Description:
% -Compute harr feature
% $Agruments
% Input;
% -clfparams: classifier parameters
% -clfparams.width: width of search window 窗口的宽
% -clfparams.height:height of search window 窗口的高
% -ftrparams: feature parameters
% -ftrparams.minNumRect: minimal number of feature rectangles
% -ftrparams.maxNumRect: maximal ....
% -M: totle number of features 特征个数
% Output:
% -px: x coordinate, size: M x ftrparms.maxNumRect
% -py: y ...
% -pw: corresponding width,size:...
% -ph: corresponding height,size:...
% -pwt:corresponding weight,size:....Range:[-1 1]
% $ History $
% - Created by Kaihua Zhang, on April 22th, 2011
%
width = clfparams.width;
height = clfparams.height;
px = zeros(M,ftrparams.maxNumRect);
py = zeros(M,ftrparams.maxNumRect);
pw = zeros(M,ftrparams.maxNumRect);
ph = zeros(M,ftrparams.maxNumRect);
pwt= zeros(M,ftrparams.maxNumRect);
%函数randi产生一个均匀分布的整数
for i=1:M
numrects = ftrparams.minNumRect + randi(ftrparams.maxNumRect-ftrparams.minNumRect)-1;
for j = 1:numrects
px(i,j) = randi(width-3); %x坐标
py(i,j) = randi(height-3); %y坐标
pw(i,j) = randi(width-px(i,j)-2); %宽
ph(i,j) = randi(height-py(i,j)-2); %高
pwt(i,j)= (-1)^(randi(2)); %权重
pwt(i,j)=pwt(i,j)/sqrt(numrects);
end
end
sampleImg.m文件
function samples = sampleImg(img,initstate,inrad,outrad,maxnum)
% $Description:
% -Compute the coordinate of sample image templates
% $Agruments
% Input;
% -img: inpute image
% -initistate: [x y width height] object position
% -inrad: outside radius of region
% -outrad: inside radius of region
% -maxnum: maximal number of samples
% Output:
% -samples.sx: x coordinate vector,[x1 x2 ...xn]
% -samples.sy: y ...
% -samples.sw: width ...
% -samples.sh: height...
% $ History $
% - Created by Kaihua Zhang, on April 22th, 2011
% - Revised by Kaihua Zhang, on May 25th, 2011
% rand('state',0);%important
inrad = ceil(inrad);
outrad= ceil(outrad);
[row,col] = size(img);
x = initstate(1);
y = initstate(2);
w = initstate(3);
h = initstate(4);
rowsz = row - h - 1;
colsz = col - w - 1;
inradsq = inrad^2;
outradsq = outrad^2;
minrow = max(1, y - inrad+1);
maxrow = min(rowsz-1, y+inrad);
mincol = max(1, x-inrad+1);
maxcol = min(colsz-1, x+inrad);
prob = maxnum/((maxrow-minrow+1)*(maxcol-mincol+1));
i = 1;
%--------------------------------------------------
%--------------------------------------------------
[r,c] = meshgrid(minrow:maxrow,mincol:maxcol);
dist = (y-r).^2+(x-c).^2;
rd = rand(size(r));
ind = (rd<prob)&(dist<inradsq)&(dist>=outradsq);
c = c(ind==1);
r = r(ind==1);
samples.sx = c';
samples.sy = r';
samples.sw = w*ones(1,length(r(:)));
samples.sh = h*ones(1,length(r(:)));
%--------------------------------------------------
% for r = minrow:maxrow
% for c = mincol:maxcol
% dist = (y-r)^2 + (x-c)^2;
% if (rand<prob)&(dist<inradsq)&(dist>=outradsq)
% samples.sx(i) = c;
% samples.sy(i) = r;
% samples.sw(i) = w;
% samples.sh(i) = h;
% i=i+1;
% end
% end
% end
integral.cpp文件，积分图像
#include <math.h>
#include "mex.h"
// compute integral img
// s(i,j) = s(i-1,j)+i(i,j)
// ii(i,j) = ii(i,j-1)+s(i,j)
// s(i,j) = s(i+j*M);
// s(0,j) = i(0,j);ii(i,0)=s(i,0)
/* Input Arguments */
#define img_IN prhs[0]
/* Output Arguments */
#define ii_OUT plhs[0]
static void integral(
double ii[],
double *img,
int M,
int N)
{
int i;
int j;
double *s = new double[M*N];
for(j=0; j<N; j++)
{
s[j*M] = img[j*M];
for(i=1; i<M; i++)
{
s[i+j*M] = s[i-1+j*M] + img[i+j*M];
}
}
for(i=0; i<M; i++)
{
ii[i] = s[i];
for(j=1; j<N; j++)
{
ii[i+j*M] = ii[i+(j-1)*M] + s[i+j*M];
}
}
delete []s;
return;
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray*prhs[] )
{
double *ii;
double *img;
mwSize M,N;
/* Check the dimensions of Y. Y can be 4 X 1 or 1 X 4. */
M = mxGetM(img_IN);
N = mxGetN(img_IN);
/* Create a matrix for the return argument */
ii_OUT = mxCreateDoubleMatrix(M, N, mxREAL);
/* Assign pointers to the various parameters */
ii = mxGetPr(ii_OUT);
img = mxGetPr(img_IN);
/* Do the actual computations in a subroutine */
integral(ii,img, M, N);
return;
}getFtrVal.m函数
function samplesFtrVal = getFtrVal(iH,samples,ftr)
% $Description:
% -Compute the features of samples
% $Agruments
% Input;
% -iH: inpute integral image
% -samples: sample templates.samples.sx:x coordinate vector, samples.sy:
% y coordinate vector
% -ftr: feature template. ftr.px,ftr.py,ftr.pw,ftr.ph,ftr.pwt
% Output:
% -samplesFtrVal: size: M x N, where M is the number of features, N is
% the number of samples
% $ History $
% - Created by Kaihua Zhang, on April 22th, 2011
% - Revised by Kaihua Zhang, 10/12/2011
sx = samples.sx;
sy = samples.sy;
px = ftr.px;
py = ftr.py;
pw = ftr.pw;
ph = ftr.ph;
pwt= ftr.pwt;
%调用函数FtrVal.mexw32，提取特征
samplesFtrVal = FtrVal(iH,sx,sy,px,py,pw,ph,pwt); %feature without preprocessing
FtrVal.cpp函数
#include <math.h>
#include "mex.h"
// compute integral img
// s(i,j) = s(i-1,j)+i(i,j)
// ii(i,j) = ii(i,j-1)+s(i,j)
// s(i,j) = s(i+j*M);
// s(0,j) = i(0,j);ii(i,0)=s(i,0)
/* Input Arguments */
/* Output Arguments */
static void getFtrVal(double samplesFtrVal[],double*iH,double *sx,double * sy,double *px, double *py, double *pw,
double *ph, double *pwt, int len_F,int len_S,int len_R,int M,int N)
{
int i,j,minJ,maxJ,minI,maxI;
int m,k;
int x,y;
int *temp = new int[len_F];
for(i=0;i<len_F; i++)
{
m=0;
for(j=0;j<len_R;j++)
{
if(px[i+j*len_F]!=0)
{
m = m+1;
}
else
{
break;
}
}
temp[i] = m;
}
for(i=0;i<len_F;i++)
for(j=0;j<len_S;j++)
{
m = 0;
x = sx[j];
y = sy[j];
for(k=0;k<temp[i];k++)
{
minJ = x-1+px[i+k*len_F];
maxJ = x-1+px[i+k*len_F]+pw[i+k*len_F]-1;
minI = y-1+py[i+k*len_F];
maxI = y-1+py[i+k*len_F]+ph[i+k*len_F]-1;
m = m+pwt[i+k*len_F]*(iH[minI+minJ*M]+iH[maxI+maxJ*M]
-iH[maxI+minJ*M]-iH[minI+maxJ*M]);
}
samplesFtrVal[i+j*len_F]=m;
}
delete []temp;
}
void mexFunction( int nlhs, mxArray *plhs[],
int nrhs, const mxArray*prhs[] )
{
double *iH,*sx,*sy,*px,*py,*pw,*ph,*pwt;
double *samplesFtrVal;
mwSize len_F,len_S,len_R,M,N;
int i,j,k,m1,n1,c;
iH = mxGetPr(prhs[0]);
sx = mxGetPr(prhs[1]);
sy = mxGetPr(prhs[2]);
px = mxGetPr(prhs[3]);//s.rect.x
py = mxGetPr(prhs[4]);//s.rect.y
pw = mxGetPr(prhs[5]);//s.rect.width
ph = mxGetPr(prhs[6]);//s.rect.height
pwt = mxGetPr(prhs[7]);//s.weight
len_F = mxGetM(prhs[3]);
len_S = mxGetN(prhs[1]);
len_R = mxGetN(prhs[3]);
M = mxGetM(prhs[0]);
N = mxGetN(prhs[0]);
plhs[0] = mxCreateDoubleMatrix(len_F,len_S, mxREAL);
samplesFtrVal = mxGetPr(plhs[0]);
getFtrVal(samplesFtrVal,iH,sx,sy,px,py,pw,ph,pwt,len_F,len_S,len_R,M,N);
return;
}classiferUpdate.m函数
function [mu1,sig1,mu0,sig0] = classiferUpdate(posx,negx,mu1,sig1,mu0,sig0,lRate)
% $Description:
% -Update the mean and variance of the gaussian classifier
% $Agruments
% Input;
% -posx: positive sample set. We utilize the posx.feature
% -negx: negative .... ... negx.feature
% -mu1: mean of positive.feature M x 1 vector
% -sig1:standard deviation of positive.feature M x 1 vector
% -mu0 : ... negative
% -sig0: ... negative
% -lRate: constant rate
% Output:
% -mu1: updated mean of positive.feature
% -sig1:... standard deviation ....
% -mu0: updated mean of negative.feature
% -sig0:.... standard variance ...
% $ History $
% - Created by Kaihua Zhang, on April 22th, 2011
% - Changed by Kaihua Zhang, on May 18th, 2011
%--------------------------------------------------
[prow,pcol] = size(posx.feature);
pmu = mean(posx.feature,2);
posmu = repmat(pmu,1,pcol);
sigm1 = mean((posx.feature-posmu).^2,2);
nmu = mean(negx.feature,2);
[nrow,ncol] = size(negx.feature);
negmu = repmat(nmu,1,ncol);
sigm0 = mean((negx.feature-negmu).^2,2);
%------------------------------------------Online MIL update method
% sig1= lRate*sig1+ (1-lRate)*sqrt(sigm1);
% mu1 = lRate*mu1 + (1-lRate)*pmu;
%
% sig0= lRate*sig0+ (1-lRate)*sqrt(sigm0);
% mu0 = lRate*mu0 + (1-lRate)*nmu;
%---------------------------------------------
%------------------------------------------Our update method
sig1= sqrt(lRate*sig1.^2+ (1-lRate)*sigm1+lRate*(1-lRate)*(mu1-pmu).^2);
mu1 = lRate*mu1 + (1-lRate)*pmu;
sig0= sqrt(lRate*sig0.^2+ (1-lRate)*sigm0+lRate*(1-lRate)*(mu0-nmu).^2);
mu0 = lRate*mu0 + (1-lRate)*nmu;ratioClassifier.m函数
function r = ratioClassifier(posx,negx,samples,M)
% $Description:
% -Compute the ratio classifier
% $Agruments
% Input;
% -posx: trained positive sample set. We utilize the posx.mu,posx.sig
% -negx: trained negative .... ... negx.mu,negx.sig
% -samples: tested samples. We utilize samples.feature
% Output:
% -r: the computed classifier with respect to samples.feature
% $ History $
% - Created by Kaihua Zhang, on April 22th, 2011
[row,col] = size(samples.feature);
mu1 = posx.mu;
sig1= posx.sig;
mu0 = negx.mu;
sig0= negx.sig;
mu1 = repmat(mu1,1,col);
sig1 = repmat(sig1,1,col);
mu0 = repmat(mu0,1,col);
sig0 = repmat(sig0,1,col);
n0= 1./(sig0+eps);
n1= 1./(sig1+eps);
e1= -1./(2*sig1.^2+eps);
e0= -1./(2*sig0.^2+eps);
x = samples.feature;
p0 = exp((x-mu0).^2.*e0).*n0;
p1 = exp((x-mu1).^2.*e1).*n1;
r = (log(eps+p1)-log(eps+p0));
数据在作者的主要上有哦。
扫一扫
3286
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
