【目标跟踪】|MOSSE原理及对应代码解释 matlab C
1原理
https://www.bilibili.com/video/av74302620/?spm_id_from=333.788.videocard.0
https://blog.csdn.net/fzp95/article/details/78385795?utm_medium=distribute.pc_relevant.none-task-blog-baidujs_title-4&spm=1001.2101.3001.4242
相關(guān)和卷積操作》
 https://blog.csdn.net/qq_17783559/article/details/82254996
1.1 相關(guān)性
 其中 f?表示 f的復(fù)共軛。correlation的直觀解釋就是衡量兩個(gè)函數(shù)在某個(gè)時(shí)刻相似程度。
1.2 MOSSE 濾波器H公式推導(dǎo)
作者提出的濾波器叫做Minimum Output Sum of Squared Error filter(MOSSE)(誤差最小平方和濾波器)。
MOSSE是一種從較少的訓(xùn)練圖像中生成濾波器的算法。
1.2.1 結(jié)論
先放結(jié)論:
 整個(gè)濾波器的模型公式:
 分子是輸入和期望輸出之間的相關(guān)性,分母是輸入的能譜。
 其中g(shù)i表示響應(yīng)輸出,fi表示輸入圖像,對(duì)多組輸入和對(duì)應(yīng)的輸出求均值,得到h,h表示濾波模板。
Q1: 對(duì)于框定的一個(gè)目標(biāo)框,怎么得到多個(gè)訓(xùn)練樣本輸入?
 A1: 可以通過任意旋轉(zhuǎn)裁剪
Q2: 訓(xùn)練時(shí)候 對(duì)應(yīng)的響應(yīng)是什么?
 A2: 一般來說,gi 可以是任何形狀。在這種情況下,gi是由ground truth生成的,它在訓(xùn)練圖像fi中有一個(gè)以目標(biāo)為中心的緊湊的(σ = 2.0)二維高斯峰.
Q3:經(jīng)過樣本訓(xùn)練得到了H ,之后怎么用?
 A3:利用得到相關(guān)性最大的,即響應(yīng)最大的地方,作為下一次跟蹤的目標(biāo)位置。
ps
在實(shí)際跟蹤的過程中我們要考慮到目標(biāo)的外觀變換等因素的影響,所以需要同時(shí)考慮目標(biāo)的m個(gè)圖像作為參考,對(duì)跟蹤框(groundtruth)進(jìn)行隨機(jī)仿射變換,獲取一系列的訓(xùn)練樣本 fi,
而 gi則是由高斯函數(shù)產(chǎn)生,并且其峰值位置是在 fi的中心位置。
獲得了一系列的訓(xùn)練樣本和結(jié)果之后,就可以計(jì)算濾波器h的值。注意:這里的f,g,h的size大小都相同。從而提高濾波器模板的魯棒性。
同時(shí),作者為了讓濾波器對(duì)與形變、光照等外界影響具有更好的魯棒性,采取了如下的模板更新策略,作者將濾波器的模型公式分為分子和分母兩個(gè)部分,每個(gè)部分都分別的進(jìn)行更新:。
 
 其中 At和 At?1分別表示的是當(dāng)前幀和上一幀的分子。
1.2.2 推導(dǎo)(可以不管)
我們需要找到一個(gè)濾波器,使其在目標(biāo)上的響應(yīng)最大,則如下公式:
 
 其中g(shù)表示響應(yīng)輸出,f表示輸入圖像,h表示濾波模板。
顯然,我們要是想獲得比較獲得響應(yīng)輸出,只需確定濾波器模板h即可。上式的計(jì)算是進(jìn)行卷積計(jì)算,這在計(jì)算機(jī)中的計(jì)算消耗時(shí)很大的,因此作者對(duì)上式進(jìn)行快速傅里葉變換(FFT)這樣卷積操作經(jīng)過FFT后就變成了點(diǎn)乘操作,極大的減少了計(jì)算量。上式變成如下形式:
 為了方便描述,將上式寫成如下形式:
后面跟蹤的任務(wù)就是找到 H?了
因?yàn)樯鲜降牟僮鞫际窃丶?jí)別的,因此要想找到,只要使其中的每個(gè)元素的MOSSE都最小即可。因此上式可轉(zhuǎn)換為如下形式:
 (w和v是H中每個(gè)元素的索引)
要想得到最小的 H*wv,只需要對(duì)上式求偏導(dǎo),并使偏導(dǎo)為0即可。即
 
 即每個(gè)位置的值為:
 
結(jié)論
整個(gè)濾波器的模型公式:
 分子是輸入和期望輸出之間的相關(guān)性,分母是輸入的能譜。
ps
但是在實(shí)際跟蹤的過程中我們要考慮到目標(biāo)的外觀變換等因素的影響,所以需要同時(shí)考慮目標(biāo)的m個(gè)圖像作為參考,對(duì)跟蹤框(groundtruth)進(jìn)行隨機(jī)仿射變換,獲取一系列的訓(xùn)練樣本 fi,
而 gi則是由高斯函數(shù)產(chǎn)生,并且其峰值位置是在 fi的中心位置。
獲得了一系列的訓(xùn)練樣本和結(jié)果之后,就可以計(jì)算濾波器h的值。注意:這里的f,g,h的size大小都相同。從而提高濾波器模板的魯棒性。
同時(shí),作者為了讓濾波器對(duì)與形變、光照等外界影響具有更好的魯棒性,采取了如下的模板更新策略,作者將濾波器的模型公式分為分子和分母兩個(gè)部分,每個(gè)部分都分別的進(jìn)行更新:。
 
 其中 At和 At?1分別表示的是當(dāng)前幀和上一幀的分子。
2 具體計(jì)算過程
首先,它需要一組訓(xùn)練圖像fiand訓(xùn)練輸出gi。
2.1 輸入目標(biāo)區(qū)域
rect
%% 選擇目標(biāo)區(qū)域 % select target from first frame im = imread(img_files(1,:)); f = figure('Name', 'Select object to track'); imshow(im); rect = getrect; close(f); clear f; center = [rect(2)+rect(4)/2 rect(1)+rect(3)/2];2.2 獲得高斯核
%% 高斯核 % plot gaussian sigma = 100; gsize = size(im);%獲取圖像尺寸 [R,C] = ndgrid(1:gsize(1), 1:gsize(2)); %兩個(gè)矩陣R和C,都是m行n列 g = gaussC(R,C, sigma, center);%通過R和C產(chǎn)生size等同于im的高斯濾波函數(shù) g = mat2gray(g); g = imcrop(g, rect); G = fft2(g); % 得到目標(biāo)區(qū)域的高斯濾波核的傅里葉變換2.3 計(jì)算 輸入的fi
其中fi是框選的目標(biāo)區(qū)域
img = imcrop(img, rect);人為地連接圖像fi的邊界會(huì)引入一個(gè)影響相關(guān)輸出的偽影。
解決方法:
首先,使用日志函數(shù)對(duì)像素值進(jìn)行轉(zhuǎn)換,這有助于低對(duì)比度的照明情況。像素值被歸一化為平均值0.0和平均值1.0。
最后,將圖像乘以一個(gè)余弦窗口,使邊緣附近的像素值逐漸減少到零。這也有一個(gè)好處,它把重點(diǎn)放在了靠近目標(biāo)中心的地方。
 原因:由于傅里葉變換是周期性的,它不尊重圖像的邊界。非周期圖像的相對(duì)邊之間的巨大不連續(xù)將導(dǎo)致有噪聲的傅里葉表示。
preprocess.m處理函數(shù)
function img = preprocess(img) [r,c] = size(img); win = window2(r,c,@hann); save win eps = 1e-5; img = log(double(img)+1); img = (img-mean(img(:)))/(std(img(:))+eps); img = img.*win; end2.4 計(jì)算初始的Hi(Ai,Bi)及生成多個(gè)訓(xùn)練模板
Ai = (G.*conj(fft2(fi)));%相關(guān)性 Bi = (fft2(fi).*conj(fft2(fi)));%fi能量譜N = 128; for i = 1:Nfi = preprocess(rand_warp(img));Ai = Ai + (G.*conj(fft2(fi)));Bi = Bi + (fft2(fi).*conj(fft2(fi))); end2.5 利用初始Hi進(jìn)行預(yù)測響應(yīng)位置,并在線更新
%預(yù)測位置Hi = Ai./Bi;fi = imcrop(img, rect); fi = preprocess(imresize(fi, [height width])); gi = uint8(255*mat2gray(ifft2(Hi.*fft2(fi))));maxval = max(gi(:))[P, Q] = find(gi == maxval);dx = mean(P)-height/2;dy = mean(Q)-width/2;rect = [rect(1)+dy rect(2)+dx width height];%更新濾波器fi = imcrop(img, rect); fi = preprocess(imresize(fi, [height width]));Ai = eta.*(G.*conj(fft2(fi))) + (1-eta).*Ai;Bi = eta.*(fft2(fi).*conj(fft2(fi))) + (1-eta).*Bi;code
matlab
%get images from source directory %% 獲取圖像 %此處僅僅用于得到圖片序列所在地址 datadir = 'H:/code/1doctor_code/matlab_code/tracking/mosse-tracker-master/data/'; dataset = 'Surfer'; path = [datadir dataset]; img_path = [path '/img/']; D = dir([img_path, '*.jpg']);%在img_path下的所有jpg后綴文件的地址放入D中seq_len = length(D(not([D.isdir])));%得到圖片總數(shù) if exist([img_path num2str(1, '%04i.jpg')], 'file')img_files = num2str((1:seq_len)', [img_path '%04i.jpg']); elseerror('No image files found in the directory.'); end %% 選擇目標(biāo)區(qū)域 % select target from first frame im = imread(img_files(1,:)); f = figure('Name', 'Select object to track'); imshow(im); rect = getrect; close(f); clear f; center = [rect(2)+rect(4)/2 rect(1)+rect(3)/2]; %% 高斯核 % plot gaussian sigma = 100; gsize = size(im);%獲取圖像尺寸 [R,C] = ndgrid(1:gsize(1), 1:gsize(2)); %兩個(gè)矩陣R和C,都是m行n列 g = gaussC(R,C, sigma, center);%通過R和C產(chǎn)生size等同于im的高斯濾波函數(shù) g = mat2gray(g); g = imcrop(g, rect); G = fft2(g); % 得到目標(biāo)區(qū)域的高斯濾波核的傅里葉變換 %% 隨機(jī)裁剪原始目標(biāo)區(qū)域,得到N個(gè)訓(xùn)練樣本 %randomly warp original image to create training set if (size(im,3) == 3) img = rgb2gray(im); end img = imcrop(img, rect);%將高斯濾波函數(shù)變換到頻域 height = size(g,1); width = size(g,2); %imresize(img, [height width])將圖片調(diào)整成濾波器的大小 fi = preprocess(imresize(img, [height width])); Ai = (G.*conj(fft2(fi))); Bi = (fft2(fi).*conj(fft2(fi)));N = 128; for i = 1:Nfi = preprocess(rand_warp(img));Ai = Ai + (G.*conj(fft2(fi)));Bi = Bi + (fft2(fi).*conj(fft2(fi))); end% MOSSE online training regimen eta = 0.25; fig = figure('Name', 'MOSSE'); t = figure; mkdir(['results_' dataset]); for i = 1:size(img_files, 1)img = imread(img_files(i,:));im = img;if (size(img,3) == 3)img = rgb2gray(img);%灰度圖endif (i == 1)Ai = eta.*Ai;Bi = eta.*Bi;%似乎沒有意義,在i=2的時(shí)候同時(shí)約了etaelse%預(yù)測位置Hi = Ai./Bi;fi = imcrop(img, rect); fi = preprocess(imresize(fi, [height width])); gi = uint8(255*mat2gray(ifft2(Hi.*fft2(fi))));maxval = max(gi(:))[P, Q] = find(gi == maxval);dx = mean(P)-height/2;dy = mean(Q)-width/2;rect = [rect(1)+dy rect(2)+dx width height];%更新濾波器fi = imcrop(img, rect); fi = preprocess(imresize(fi, [height width]));Ai = eta.*(G.*conj(fft2(fi))) + (1-eta).*Ai;Bi = eta.*(fft2(fi).*conj(fft2(fi))) + (1-eta).*Bi;end% visualizationtext_str = ['Frame: ' num2str(i)];box_color = 'green';position=[1 1];result = insertText(im, position,text_str,'FontSize',15,'BoxColor',...box_color,'BoxOpacity',0.4,'TextColor','white');result = insertShape(result, 'Rectangle', rect, 'LineWidth', 3);imwrite(result, ['results_' dataset num2str(i, '/%04i.jpg')]);imshow(result);drawnow;rect endc
// This file ispart of the OpenCV project. // It is subject to the license terms in the LICENSEfile found in the top-level directory // of this distribution and athttp://opencv.org/license.html. // //[1] David S. Bolme et al. "Visual Object Trackingusing Adaptive Correlation Filters" // http://www.cs.colostate.edu/~draper/papers/bolme_cvpr10.pdf // // // credits: // Kun-Hsin Chen: for initial c++ code // Cracki: for the idea of only converting the usedpatch to gray // #include "opencv2/tracking.hpp" namespace cv { namespace tracking { struct DummyModel :TrackerModel { virtual void modelUpdateImpl(){} virtual void modelEstimationImpl( const std::vector<Mat>& ){} }; const double eps=0.00001; // fornormalization const double rate=0.2; //learning rate const double psrThreshold=5.7; //no detection, if PSR is smaller than this struct MosseImpl :TrackerMOSSE { protected: Point2d center;//center of the bounding box Size size; //size ofthe bounding box Mat hanWin; Mat G; //goal Mat H, A,B; //state// Element-wisedivision of complex numbers in src1 and src2 Mat divDFTs( const Mat &src1, const Mat &src2 ) const { Mat c1[2],c2[2],a1,a2,s1,s2,denom,re,im; // split into re and im per src cv::split(src1, c1); cv::split(src2, c2); // (Re2*Re2 + Im2*Im2) = denom // denom is same forboth channels cv::multiply(c2[0], c2[0], s1); cv::multiply(c2[1], c2[1], s2); cv::add(s1,s2, denom); // (Re1*Re2 + Im1*Im1)/(Re2*Re2 + Im2*Im2) = Re cv::multiply(c1[0], c2[0], a1); cv::multiply(c1[1], c2[1], a2); cv::divide(a1+a2, denom, re, 1.0 ); // (Im1*Re2 - Re1*Im2)/(Re2*Re2 + Im2*Im2) = Im cv::multiply(c1[1], c2[0], a1); cv::multiply(c1[0], c2[1], a2); cv::divide(a1+a2, denom, im, -1.0); // Merge Re and Im back into a complex matrix Mat dst,chn[] = {re,im}; cv::merge(chn, 2, dst); return dst; } void preProcess( Mat &window) const { window.convertTo(window, CV_32F); log(window+ 1.0f, window); //normalize Scalarmean,StdDev; meanStdDev(window, mean, StdDev); window =(window-mean[0]) /(StdDev[0]+eps); //Gaussain weighting window =window.mul(hanWin); } double correlate( const Mat &image_sub, Point&delta_xy ) const//計(jì)算相對(duì)位移 { MatIMAGE_SUB, RESPONSE, response; // filter in dft space dft(image_sub, IMAGE_SUB, DFT_COMPLEX_OUTPUT); mulSpectrums(IMAGE_SUB, H, RESPONSE, 0, true ); idft(RESPONSE, response, DFT_SCALE|DFT_REAL_OUTPUT); // update center position double maxVal; Point maxLoc; minMaxLoc(response, 0, &maxVal, 0, &maxLoc); delta_xy.x= maxLoc.x - int(response.size().width/2); delta_xy.y= maxLoc.y - int(response.size().height/2); // normalize response Scalarmean,std; meanStdDev(response, mean, std); return (maxVal-mean[0]) / (std[0]+eps); // PSR } Mat randWarp( const Mat& a ) const { cv::RNGrng(8031965); // random rotation double C=0.1; double ang = rng.uniform(-C,C); double c=cos(ang), s=sin(ang); // affine warp matrix Mat_<float> W(2,3); W <<c + rng.uniform(-C,C), -s + rng.uniform(-C,C), 0, s +rng.uniform(-C,C), c +rng.uniform(-C,C), 0;// random translation Mat_<float> center_warp(2, 1); center_warp << a.cols/2, a.rows/2; W.col(2) = center_warp - (W.colRange(0, 2))*center_warp; Mat warped;warpAffine(a, warped, W, a.size(), BORDER_REFLECT); return warped; } virtual bool initImpl( const Mat& image, const Rect2d& boundingBox ){ model =makePtr<DummyModel>(); Mat img; if (image.channels() == 1) img =image; else cvtColor(image, img, COLOR_BGR2GRAY); int w = getOptimalDFTSize(int(boundingBox.width)); int h = getOptimalDFTSize(int(boundingBox.height)); //Get the center position int x1 = int(floor((2*boundingBox.x+boundingBox.width-w)/2)); int y1 = int(floor((2*boundingBox.y+boundingBox.height-h)/2)); center.x =x1 + (w)/2; center.y =y1 + (h)/2; size.width= w; size.height= h; Mat window;getRectSubPix(img, size, center, window); createHanningWindow(hanWin, size, CV_32F); // goal Matg=Mat::zeros(size,CV_32F); g.at<float>(h/2, w/2) = 1;GaussianBlur(g, g, Size(-1,-1),2.0); double maxVal; minMaxLoc(g, 0,&maxVal); g = g /maxVal; dft(g, G,DFT_COMPLEX_OUTPUT); // initial A,B and H A =Mat::zeros(G.size(), G.type()); B =Mat::zeros(G.size(), G.type()); for(int i=0; i<8; i++) { Matwindow_warp = randWarp(window); preProcess(window_warp); MatWINDOW_WARP, A_i, B_i; dft(window_warp, WINDOW_WARP, DFT_COMPLEX_OUTPUT); mulSpectrums(G ,WINDOW_WARP, A_i, 0, true); mulSpectrums(WINDOW_WARP, WINDOW_WARP, B_i, 0, true); A+=A_i;B+=B_i;} H =divDFTs(A,B); return true; } virtual bool updateImpl( const Mat& image,Rect2d& boundingBox ) { if (H.empty()) // not initialized return false; Matimage_sub; getRectSubPix(image, size, center, image_sub); if (image_sub.channels() != 1) cvtColor(image_sub, image_sub, COLOR_BGR2GRAY); preProcess(image_sub); Pointdelta_xy; double PSR =correlate(image_sub, delta_xy); if (PSR < psrThreshold) return false; //update location center.x +=delta_xy.x; center.y +=delta_xy.y; Matimg_sub_new; getRectSubPix(image, size, center, img_sub_new); if (img_sub_new.channels() !=1) cvtColor(img_sub_new, img_sub_new, COLOR_BGR2GRAY); preProcess(img_sub_new); // new state for A and B Mat F, A_new,B_new; dft(img_sub_new, F, DFT_COMPLEX_OUTPUT); mulSpectrums(G, F, A_new, 0, true ); mulSpectrums(F, F, B_new, 0, true ); // update A ,B, and H A = A*(1-rate) + A_new*rate; B = B*(1-rate) + B_new*rate; H =divDFTs(A, B); // return tracked rect double x=center.x, y=center.y; int w = size.width,h=size.height; boundingBox= Rect2d(Point2d(x-0.5*w,y-0.5*h), Point2d(x+0.5*w, y+0.5*h)); return true; } public: MosseImpl() {isInit = 0; } // dummy implementation. virtual void read( const FileNode& ){} virtual void write( FileStorage& ) const{} }; // MosseImpl } // tracking Ptr<TrackerMOSSE> TrackerMOSSE::create() { returnmakePtr<tracking::MosseImpl>(); } } // cv總結(jié)
以上是生活随笔為你收集整理的【目标跟踪】|MOSSE原理及对应代码解释 matlab C的全部內(nèi)容,希望文章能夠幫你解決所遇到的問題。
 
                            
                        - 上一篇: 翻译记忆软件:Trados 7/2006
- 下一篇: jenkins 部署文档
