matlab練習程序（直方圖反向投影）

本文轉載自查看原文 2012-11-22 13:21 5701 Matlab練習程序

　　做meanshift物體跟蹤的時候中間有一步叫做直方圖反向投影，所以我就先實現了這樣一個步驟。

　　直方圖反向投影說白了就是模板匹配，給定一個較小的目標模板，然后再逐個遍歷原圖像和模板圖像相同的圖像塊的，對比圖像塊和模板的直方圖，然后把比較結果存入一個新的圖像中，新圖像中的全局極值就是模板在原圖像中所在的位置。這里主要麻煩的是怎么比較兩個圖像塊的直方圖，Opencv中實現了5種對比的方法，所以我在這里也對應的實現了5種方法。

　　5種方法分別是correl（相關）、chisqr（卡方）、intersect（相交）、bhattacharyya（名字很長--!!）、emd（earth mover's distance）。

　　下面是這5種方法的公式，H1是第一個圖像塊的直方圖序列，H2是第二的圖像塊的直方圖序列：

　　相關法：

　　其中：

　　卡方法：

　　相交法：

　　bhattacharyya:

　　emd:

　　前4種方法比較簡單，emd就比較麻煩了，下面是我自己的理解。這個算法使用原序列生成f和d這兩個矩陣，然后對這兩個矩陣進行處理。d矩陣按我的理解就是H1中各個元素和H2中各個元素在位置上的距離，比如H1和H2都是只有三個元素的序列，H1(1)和H2(1)的距離就是0，所以d(1,1)=0，H1(3)和H2(2)的距離是1，所以d(3,2)=1，以此d=[0 1 2;1 0 1;2 1 0]，這里d和H1，H2中的值沒關系。f的確定就比較麻煩了，看看約束條件就讓人頭疼了。按照約束條件是沒辦法寫程序的，不過我想到一個巧妙的方法求的的結果正好滿足那四個約束條件：首先不管前三個條件，只使用第四個條件求得矩陣f，然后對f從最右列開始依次減去相鄰的左邊的列，然后對f從最下行開始依次減去相鄰的上邊的行，求得的結果正好滿足條件。至於為什么這樣我還不得其解，原因就指望哪位數學大牛去探索吧。

　　下面是代碼：

　　main.m

close all;
clear all;
clc;

img=imread('lena.jpg');
imshow(img);
[m n]=size(img);

w=imcrop();      %這里把要裁剪的圖像框出來
[H W]=size(w);
hist1=histcount(w);

HH=floor(H/2);
WW=floor(W/2);

imgn=zeros(m+2*HH+1,n+2*WW+1);
imgn(HH+1:m+HH,WW+1:n+WW)=img;
imgn(1:HH,WW+1:n+WW)=img(1:HH,1:n); 
imgn(1:m+HH,n+WW+1:n+2*WW+1)=imgn(1:m+HH,n:n+WW);
imgn(m+HH+1:m+2*HH+1,WW+1:n+2*WW+1)=imgn(m:m+HH,WW+1:n+2*WW+1);
imgn(1:m+2*HH+1,1:WW)=imgn(1:m+2*HH+1,WW+1:2*WW);

re1=imgn;
re2=imgn;
re3=imgn;
re4=imgn;
re5=imgn;
for i=HH+1:m+HH
    for j=WW+1:n+WW
        s=imgn(i-HH:i+HH,j-WW:j+WW);
        hist2=histcount(s);
     
        re1(i,j)=correl(hist1,hist2);           %相關法
        re2(i,j)=chisqr(hist1,hist2);           %卡方法
        re3(i,j)=intersect(hist1,hist2);        %相交法
        re4(i,j)=bhattacharyya(hist1,hist2);    %名字很長的法
        re5(i,j)=emd(hist1,hist2);     %由於沒有優化，速度實在太慢了，至少運行一晚上，慎用！！
    end
end
figure;
re1=re1(HH+1:m+HH,WW+1:n+WW);
imshow(mat2gray(re1));

figure;
re2=re2(HH+1:m+HH,WW+1:n+WW);
imshow(mat2gray(re2));

figure;
re3=re3(HH+1:m+HH,WW+1:n+WW);
imshow(mat2gray(re3));

figure;
re4=re4(HH+1:m+HH,WW+1:n+WW);
imshow(mat2gray(re4));

figure;
re5=re5(HH+1:m+HH,WW+1:n+WW);
imshow(mat2gray(re5));

　　histcount.m 統計直方圖

function hist=histcount(w)
    [H W]=size(w);
    w=uint8(w);
    hist=zeros(1,256);
    for i=1:H
        for j=1:W
            hist(w(i,j)+1)=hist(w(i,j)+1)+1;
        end
    end

    hist=hist/(H*W);

end

　　correl.m 相關法

function d=correl(H1,H2)
    d=sum(H1-mean(H1).*(H2-mean(H2)))/sqrt(sum((H1-mean(H1)).^2)*sum((H2-mean(H2)).^2));
end

　　chisqr.m 卡方法

function d=chisqr(H1,H2)
    d=sum(((H1-H2).^2)/(H1+H2));
end

　　intersect.m 相交法

function d=intersect(H1,H2)    
    d=sum(min(H1,H2));
end

　　bhattacharyya.m

function d=bhattacharyya(H1,H2)
    d=sqrt(1-sum(sqrt(H1.*H2))/sqrt(sum(H1)*sum(H2)));
end

　　emd.m 這里沒優化，慎用！

function re=emd(H1,H2)
    m=length(H1);
    n=length(H2);
    f=zeros(m,n);
    d=zeros(m,n);
    
    for i=1:m
        for j=1:n
            if i==j
                d(i,j)=0;
            end       
            if j>i
                d(i,j)=j-i;
            end       
            if j<i
                d(i,j)=i-j;
            end    
            f(i,j)=min(sum(H1(1:i)),sum(H2(1:j)));              
        end   
    end
    
    for i=m:-1:2
        f(:,i)=f(:,i)-f(:,i-1);
    end
    
    for j=n:-1:2
        f(j,:)=f(j,:)-f(j-1,:);
    end

    re=(sum(sum(f.*d)))/sum(sum(f));
    
end