压缩感知跟踪（一）

  闲着也是闲着，不如看看代码吧，最近开始接触压缩感知跟踪，是从朋友那里看到的，所以想找了一些代码跑跑看看，算是体验一下吧！！

  论文：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: [e-mail]% Date: 11/12/2011% Revised by Kaihua Zhang, 15/12/2011% Revised by Kaihua Zhang, 11/7/2012clc;clear all;close all;%----------------------------------rand('state',0);%----------------------------------%添加路径addpath('./data');%----------------------------------%加载文件load init.txt;%[x y width height]initstate = init;%initial tracker%----------------------------Set pathimg_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~8trparams.init_postrainrad = 4.0;%radical scope of positive samples%目标初始状态[x,y,w,h]trparams.initstate = initstate;% object position [x y width height]%新的一帧中的搜索窗口的半径，通常设置为15~35trparams.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-6ftrparams.maxNumRect = 4;%-------------------------%弱分类器个数M = 50;% number of all weaker classifiers, i.e,feature pool%-------------------------%均值posx.mu = zeros(M,1);% mean of positive featuresnegx.mu = zeros(M,1);%方差posx.sig= ones(M,1);% variance of positive featuresnegx.sig= ones(M,1);%-------------------------Learning rate parameter%学习速率，通常设置为0.7~0.95lRate = 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.mexw32iH = integral(img);%Compute integral imageposx.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 cornery = initstate(2);w = initstate(3);% width of the rectangleh = 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 parametersend

HaarFtr.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);%importantinrad = 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 */#defineimg_INprhs[0]/* Output Arguments */#defineii_OUTplhs[0]static void integral(   doubleii[],   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/2011sx = 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.xpy = mxGetPr(prhs[4]);//s.rect.ypw = mxGetPr(prhs[5]);//s.rect.widthph = mxGetPr(prhs[6]);//s.rect.heightpwt = mxGetPr(prhs[7]);//s.weightlen_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 methodsig1= 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));

数据在作者的主要上有哦。