關於DPM(Deformable Part Model)算法中模型可視化的解釋

搭建了自己的博客平台，本文地址：http://masikkk.com/blog/DPM-model-visualization/

        DPM源代碼(voc-release)中的模型可視化做的還算相當炫酷的，能夠讓我們直觀的看到訓練好的模型，甚至我們不用去做模型的評價，直接依據肉眼的觀察。就能大致了解一個目標訓練的好不好，比方我訓練一個人體模型。那他的可視化圖當然就是越接近人體越好。

        以下是對DPM源代碼中有關模型可視化部分代碼的分析，通過分析這些代碼。有助於更好的理解DPM模型。

        注意：我的源代碼版本號是voc-release3.1，第4版往后的模型變得更復雜，這里不討論。

        有關模型可視化的代碼主要在visualizemodel.m，foldHOG.m和HOGpicture.m中。

（1）簡化分類器參數向量(或者叫濾波器權重向量)

        DPM中使用降維后的31維HOG特征向量，所以，與之相應的，訓練好的模型的參數向量也是31維的，為了方便可視化。須要將31維的參數向量簡化為9維，foldHOG函數就負責參數向量的簡化。

        31維的HOG特征向量是有分段含義的，例如以下：

        設C是聚合有9個對照度不敏感方向的像素級特征映射而獲得的基於cell的特征映射。D是聚合有18個對照度敏感方向的像素級特征而獲得的基於cell的特征映射。用4種不同的歸一化方法對C(i,j)和D(i,j)進行歸一化和截斷(限幅)。能夠獲得一個4*(9+18)=108維的特征向量F(i,j)。實際中我們使用此108維向量的一個解析投影，此投影由以下幾個統計量定義：27個在不同歸一化因子上的累加和(即列的和)。F中的每一個方向通道相應一個；以及4個在不同方向(9維對照度不敏感方向)上的累加和(即行的和)，每一個歸一化因子相應一個。cell尺寸k=8。截斷(限幅)閾值α=0.2。

終於的特征映射是31維向量G(i,j)。當中27維相應不同的方向通道(9個對照度不敏感方向和18個對照度敏感方向)，剩下4維表示(i,j)周圍4個cell組成的block的梯度能量。

        所以，foldHOG中的簡化過程就是將31維參數向量的后4維丟棄，將前27維進行負值抑制和折疊累加，縮減為9維的參數向量。

        foldHOG函數源代碼凝視例如以下：

function f = foldHOG(w)
% 簡化濾波器向量w。用以可視化顯示模型
% 將 width*height*31 的濾波器參數向量w濃縮為 width*height*9 的向量(width和height是濾波器的寬度和高度)
% 返回值f是一個width*height*9的矩陣
%
% f = foldHOG(w)
% Condense HOG features into one orientation histogram.
% Used for displaying a feature.


% max(w(:,:,1:9),0)返回w(:,:,1:9)中元素和0兩者中的較大值(去除負權重)。返回結果組成一個width*height*9維的矩陣
% 所以以下的處理相當於把濾波器參數向量w沿第三維維折疊了兩次，形成一個了一個width*height*9的簡化版的參數向量
% w的第三維的長度為31，僅僅使用了前27個值，舍棄了后面4個值。這和DPM中使用的31維HOG特征向量所代表的意義有關。
% DPM中的特征向量為31維，當中前27維相應不同的方向通道(9個對照度不敏感方向和18個對照度敏感方向)，
% 剩下4維表示(i,j)周圍4個cell組成的block的梯度能量。

f = max(w(:,:,1:9),0) + max(w(:,:,10:18),0) + max(w(:,:,19:27),0);

（2）生成濾波器權重向量的可視化圖

        這一工作在HOGpicture函數中完畢，此函數負責為簡化后的w*h*9維的權重向量生成可視化圖。

        首先要生成一個間隔為20度的方向坐標基。然后將濾波器向量中的點(i,j,k)向坐標基中的方向k上投影，而用該點的值衡量(i,j)在k方向的幅度。

        HOGpicture函數源代碼凝視例如以下：

function im = HOGpicture(w, bs)
% 畫出HOG正權重w的條紋圖像 
% 參數：
% w：簡化后的width*height*9的HOG正權重向量(width和height是濾波器的寬度和高度)
% bs：生成的圖像im相比於濾波器尺寸的擴大倍數
% 返回值：
% im：濾波器權重向量的可視化圖，是大小為(width*bs)*(height*bs)的圖像
%
% HOGpicture(w, bs)
% Make picture of positive HOG weights.

% 為間隔20度的9個方向生成條紋線
% 事實上bim相當於有9個方向的方向空間的一個坐標基。將濾波器向量中某點的值w(i,j,:)向bim的各個方向投影能夠反映每一個點的方向分布
% construct a "glyph" for each orientaion
bim1 = zeros(bs, bs); % 生成一個bs*bs的全零矩陣bim1
bim1(:,round(bs/2):round(bs/2)+1) = 1; % 將bim1的中間兩個豎條的值置為1
bim = zeros([size(bim1) 9]); % 生成一個bs*bs*9的全零矩陣bim。能夠將bim看做9層bim1疊加在一起
bim(:,:,1) = bim1; % 將bim的第1層bim(:,:,1)賦值為bim1
% 接下來通過20度遞進的順時針旋轉依次生成bim的第2到9層，並將旋轉后的矩陣裁剪為和bim1同樣大小
% 比如，i=2時，順時針旋轉bim1，並裁剪為和bim1同樣的大小，賦值給bim的第2層
for i = 2:9,
  bim(:,:,i) = imrotate(bim1, -(i-1)*20, 'crop'); % 依次順時針旋轉20度，將結果賦值給bim的第2到9層
end


% 通過加入帶有方向權重的條紋來繪制正權重的可視化圖
% make pictures of positive weights bs adding up weighted glyphs
s = size(w); % height * width * 9
w(w < 0) = 0; % 保證w中所有是正權重
im = zeros(bs*s(1), bs*s(2)); % 生成一個(width*bs)*(height*bs)的圖像，即將濾波器w的尺寸擴大bs倍

% 遍歷濾波器權重向量w，將每一個坐標的值投影到9個方向上，然后擴大bs倍畫到圖像im上
for i = 1:s(1), % 第i行(w原尺寸)
  iis = (i-1)*bs+1:i*bs; % 相應在圖像im上的橫坐標
  for j = 1:s(2), % 第j列(w原尺寸)
    jjs = (j-1)*bs+1:j*bs; % 相應在圖像im上的縱坐標
    for k = 1:9, % 遍歷9個方向
        % bim(:,:,k) * w(i,j,k)：假設濾波器向量在方向k上有正值的話，將這個值w(i,j,k)投影到大小為bs*bs的方向坐標基bim的方向k上，
        % 然后將9個方向上的投影累加，累加值反映了濾波器中(i,j)位置在各個方向上的幅度大小，最后將累加值放到擴大bs倍的顯示圖像im的相應位置上
        im(iis,jjs) = im(iis,jjs) + bim(:,:,k) * w(i,j,k); 
    end
    %imagesc(im); % 自己加入的語句。分析代碼用，顯示繪圖過程
  end
  %imagesc(im); % 自己加入的語句。調試用。顯示繪圖過程
end

以下圖1-6是HOGpicture中一個根濾波器可視化圖的繪制過程：

                     圖1。根濾波器權重向量點(1,1)的可視化                                                    圖2，加上點(1,2)的可視化

                                圖3，完畢第1行的可視化                                                                圖4。完畢前2行的可視化

                  圖5。完畢前10行的可視化。能夠看出人形了                                          圖6，完畢整個根濾波器的可視化      

（3）在visualizemodel函數中進行一些后處理。切割畫圖區域。依次調用HOGpicture畫出根濾波器和各個部件濾波器的可視化圖。以及各個部件的變形花費圖。

        在visualizemodel中調用HOGpicture畫出根濾波器的可視化圖。返回值為圖6，然后將像素值擴充到[0,255]並轉換為8位無符號整型，得到圖7：

                         圖7，根濾波器可視化圖_Uint8

        然后切割畫圖區域，將根濾波器的可視化圖畫到指定區域，如圖8；再轉換為灰度圖，如圖9。

                                圖8，根濾波器_subplot                                                                            圖9，根濾波器_subplot_gray

        再之后依次繪制各個部件的可視化圖。如圖10是人體頭部的可視化圖，覆蓋到根濾波器的相應位置。

       圖10，頭部部件可視化

         圖11，左圖是根濾波器，右圖是各個部件的可視化圖覆蓋到根的相應位置后的可視化圖

        從圖11中能夠看出，部件濾波器明顯要比根濾波器仔細，能提供很多其它細節。

        最后就是生成各個部件的變形花費圖了，這要用到各個部件的變形信息，在模型的defs[]數組中。defs數組中。每一個部件相應一個錨點坐標和一個變形花費參數(4維向量)。

計算部件內每一個位置距離部件中心的距離，用變形特征向量v = [Δx^2, Δx, Δy^2, Δy]' 和 部件的變形花費 相乘，得到的結果能夠反映此位置的變形花費。值越大。說明變形費用越高，表明不是部件的理想位置；值越小，說明變形費用越低。表明是該部件的理想位置。

反映到變形花費圖上，越亮(白)的地方花費越大。越暗(黑)的地方花費越小。

        例如以下圖12是頭部的變形花費圖：

                      圖12，頭部的變形花費圖

        最后。就獲得了完整的模型可視化圖

圖13。左圖：根濾波器的可視化圖。中圖：各個部件的可視化圖覆蓋到根的相應位置后的可視化圖；右圖：各個部件的變形花費圖

        visualizemodel函數源代碼凝視例如以下：

function visualizemodel(model, components)
% 繪制模型的可視化圖像
% 參數：
% model：要可視化的模型
% components：指定可視化某個組件模型
%
% visualizemodel(model)
% Visualize a model.

clf; % 刪除當前畫圖
if nargin < 2 % 未指定可視化哪個組件模型，則可視化全部組件模型
  components = 1:model.numcomponents;
end

% 依次可視化每一個組件模型
k = 1;
for i = components
  visualizecomponent(model, i, length(components), k);
  k = k+1;
end


% 可視化一個組件模型
% 參數
% model：要可視化的模型
% c：當前要可視化的第c個組件
% nc：此模型中總的組件個數
% k：指定繪制區域
function visualizecomponent(model, c, nc, k)

pad = 2; % 填充寬度
bs = 20; % 生成的圖像相比於濾波器尺寸的擴大倍數

% 將 width*height*31 的濾波器參數向量濃縮為 width*height*9 的向量，所以返回值w是一個 width*height*9 的向量(width和height是濾波器的寬度和高度)
w = foldHOG(model.rootfilters{model.components{c}.rootindex}.w); % 簡化組件c的根濾波器向量，用以可視化顯示模型

scale = max(w(:)); % w(:)返回由width*height*9的矩陣w的全部元素組成的一維向量，所以scale是w中全部元素的最大值，scale是一個標量

im = HOGpicture(w, bs); % 畫出濾波器權重向量w的可視化圖，返回值im是大小為(width*bs)*(height*bs)的圖像，像素值為double型
%imagesc(im); % 自己加入的語句，分析代碼用。顯示HOGpicture的畫圖結果

im = imresize(im, 2); % 將im的尺寸擴大一倍
%imagesc(im); % 自己加入的語句，分析代碼用，顯示圖像

im = padarray(im, [pad pad], 0); % 填充圖像邊界，在圖像im的上下左右各填充pad行(列)零值。填充后im的大小為(width*bs+pad*2)*(height*bs+pad*2)
im = uint8(im * (255/scale)); % 將im的值擴充到[0-255]並轉換為8位無符號整型
%imagesc(im); % 自己加入的語句，分析代碼用。顯示圖像

% 切割畫圖區域並畫出根濾波器的可視化圖
numparts = length(model.components{c}.parts); % numparts：組件c的部件個數
% 依據組件個數和是否含有部件來切割畫圖區域
% 對於含nc個組件的模型，將畫圖區域分為nc行
if numparts > 0 % 對於有部件的模型，再將每行分為3列
  subplot(nc,3,1+3*(k-1)); % 選中切割后的第1+3(k-1)個畫圖區域。即每行的第一列
else % 對於沒有部件的模型，每行僅僅有一列
  subplot(nc,1,k); % 選中每行第一列的畫圖區域
end
imagesc(im); % 縮放數據並顯示為圖片。畫到上一步選中的畫圖區域中
colormap gray; % 設為灰度圖
axis equal; % 使橫縱坐標的刻度同樣
axis off; % 不顯示坐標軸


% 畫出組件c的帶部件的可視化圖im和變形花費圖def_im
% draw parts and deformation model
if numparts > 0 % 僅僅對含有部件的模型進行以下的操作
  
  def_im = zeros(size(im)); % 初始化變形花費圖def_im，和根濾波器的可視化圖大小同樣
  def_scale = 500;
  
  % 遍歷組件c的各個部件
  for i = 1:numparts
    
    % 生成部件i的可視化圖
    w = model.partfilters{model.components{c}.parts{i}.partindex}.w; % 組件c的第i個部件的濾波器向量，尺寸為：width*height*31
    p = HOGpicture(foldHOG(w), bs); % 簡化濾波器向量w為width*height*9，並生成可視化圖，返回值p是大小為(width*bs)*(height*bs)的圖像，像素值為double型
    %clf;imagesc(p); % 自己加入的語句。分析代碼用，顯示圖像
    
    p = padarray(p, [pad pad], 0); % 填充圖像邊界。在圖像p的上下左右各添加pad行(列)零值。填充后p的大小為(width*bs+pad*2)*(height*bs+pad*2)
    p = uint8(p * (255/scale)); % 將p的值擴充到[0-255]並轉換為8位無符號整型
    %imagesc(p); % 自己加入的語句。分析代碼用，顯示圖像
    
    % 將部件i的可視化圖p的上下左右寬度為pad*2區域的邊界的值設為128，也就是加邊框
    p(:,1:2*pad) = 128;
    p(:,end-2*pad+1:end) = 128;
    p(1:2*pad,:) = 128;
    p(end-2*pad+1:end,:) = 128;
    %imagesc(p); % 自己加入的語句。分析代碼用，顯示圖像
    
    % paste into root 將部件i的可視化圖p覆蓋到根濾波器可視化圖圖像im的相應位置上
    def = model.defs{model.components{c}.parts{i}.defindex}; % 組件c的第i個部件的錨點信息
    x1 = (def.anchor(1)-1)*bs+1; % 部件i的錨點(左上角點)相應在根濾波器的可視化圖im中的坐標
    y1 = (def.anchor(2)-1)*bs+1; 
    x2 = x1 + size(p, 2)-1; % 部件i的右下角點在根濾波器的可視化圖im中的坐標
    y2 = y1 + size(p, 1)-1;
    im(y1:y2, x1:x2) = p; % 覆蓋到根濾波器的指定位置
    %imagesc(p); % 自己加入的語句，分析代碼用。顯示圖像
    
    % deformation model 生成部件i的變形花費圖並拷貝到總體的變形花費圖def_im中
    probex = size(p,2)/2; % p的寬度的一半
    probey = size(p,1)/2; % p的高度的一半
    % 生成p的每一個位置(忽略邊框)的變形花費值
    for y = 2*pad+1:size(p,1)-2*pad % 第y行
      for x = 2*pad+1:size(p,2)-2*pad % 第x列
        px = ((probex-x)/bs); % 點(y,x)距部件i中心的水平距離，Δx
        py = ((probey-y)/bs); % 點(y,x)距部件i中心的垂直距離。Δy
        v = [px^2; px; py^2; py]; % 偏移量及其平方組合成一個變形特征向量v = [Δx^2, Δx, Δy^2, Δy]'
        % 變形特征向量 乘以 變形花費參數 並進行縮放。得到點(y,x)的變形花費值，保存在p中
        % 依據計算公式可知，距離部件i中心越遠的地方，變形花費越大，反應到變形花費圖中就是越亮的地方變形花費越大
        p(y, x) = def.w * v * def_scale; 
      end
    end
    def_im(y1:y2, x1:x2) = p; % 將p中保存的部件i各個位置的變形花費值拷貝到總體的變形花費圖def_im的相應位置
  end
  
  % 在相應畫圖區域畫出組件c的帶部件的可視化圖im
  % plot parts
  subplot(nc,3,2+3*(k-1)); % 選中每行的第2列
  imagesc(im); % 繪制圖片
  colormap gray; % 設為灰度圖
  axis equal; % 使橫縱坐標的刻度同樣
  axis off; % 不顯示坐標軸
  
  % 在相應畫圖區域畫出組件c的各個部件的變形花費圖
  % plot deformation model
  subplot(nc,3,3+3*(k-1)); % 選中每行的第3列
  imagesc(def_im); % 繪制圖片
  colormap gray; % 設為灰度圖
  axis equal; % 使橫縱坐標的刻度同樣
  axis off; % 不顯示坐標軸
end

% set(gcf, 'Color', 'white')

相關鏈接：

Deformable Part Model 相關網頁：http://www.cs.berkeley.edu/~rbg/latent/index.html

Pedro Felzenszwalb的個人主頁：http://cs.brown.edu/~pff/

PASCAL VOC 目標檢測挑戰：http://pascallin.ecs.soton.ac.uk/challenges/VOC/

A Discriminatively Trained, Multiscale, Deformable Part Model [CVPR 2008] 中文翻譯

Object Detection with Discriminatively Trained Part Based Models [PAMI 2010]中文翻譯

有關可變形部件模型(Deformable Part Model)的一些說明

在Windows下執行Felzenszwalb的Deformable Part Models(voc-release4.01)目標檢測matlab源代碼

在Windows下執行Felzenszwalb的star-cascade DPM(Deformable Part Models)目標檢測Matlab源代碼

在windows下執行Felzenszwalb的Deformable Part Model(DPM)源代碼voc-release3.1來訓練自己的模型

用DPM(Deformable Part Model，voc-release3.1)算法在INRIA數據集上訓練自己的人體檢測模型