DPM中的HOG源码的Matlab版重写-《小超教你写论文》系列第四部分第一篇

DPM中的HOG源码的Matlab版重写-《小超教你写论文》系列第四部分第一篇

《小超教你写论文》系列前三部分分别翻译了一篇文章；对文章中的公式进行了推导；介绍了作者使用的数据库。作为系列的第四部分，开始对原文中的想法进行实现。

前一段陷入了一点误区，一直在思考一个算法，总想把一切都弄得非常明白，耗费了很多时间。其实后来发现没有必要，我们的目的是写文章，而写文章需要我们实现自己的算法。如果我们的算法有从别人那里借鉴的部分，只需要将相同的部分弄明白就可以了，没必要完整彻底地理解别人的算法（ps：这里是指以写文章为目的的过程，并非所有科研过程）。

算法中，用到了DPM中的HOG特征，所以需要实现它。但在voc-release3.1的源码中，我发现HOG的实现是C语言版的，原大神可能是为了加快运算速度，但自己用起来不是很方便，所以就按照C语言写了一个Matlab版，几乎完全一样。其实，其中梯度方向分区的部分我曾经写过速度更快的方法，但为了与原版保持一致，没有做改变。

具体代码如下，包含了HOG特征求取和HOG特征图像化两个部分。

function feat = HOG(input,sbin)
tic
uu = [1.0000;
0.9397;
0.7660;
0.500;
0.1736;
-0.1736;
-0.5000;
-0.7660;
-0.9397];

vv = [0.0000;
0.3420;
0.6428;
0.8660;
0.9848;
0.9848;
0.8660;
0.6428;
0.3420];

%读入图片，若为灰度，将灰度转化为彩色
input = imread('6.jpg');
%figure,imshow(input);
sbin = 8;
%input = color(input);
im = input;
blocks = zeros(2,1);
dims = zeros(3,1);
[dims(1),dims(2),dims(3)] = size(input);

blocks(1) = round(dims(1)/sbin);%行
blocks(2) = round(dims(2)/sbin);%列
hist = zeros(blocks(1),blocks(2),18);
norm = zeros(blocks(1),blocks(2));

out = zeros(3,1);
out(1) = max(blocks(1)-2,0);%减2是为了方便一会儿每个block利用旁边四个block中进行归一化
out(2) = max(blocks(2)-2,0);
out(3) = 27+4;
feat = zeros(out(1),out(2),out(3));%用来储存最后特征

visible = zeros(2,1);
visible(1) = blocks(1)*sbin;
visible(2) = blocks(2)*sbin;
%既然原代码是先列再行，我们也先列再行
for x = 2:visible(2)-1
for y = 2:visible(1)-1
%第一通道
dy = im(y+1,x,1)-im(y-1,x,1);
dx = im(y,x+1,1)-im(y,x-1,1);
v = dy*dy+dx*dx;
%第二通道
dy2 = im(y+1,x,2)-im(y-1,x,2);
dx2 = im(y,x+1,2)-im(y,x-1,2);
v2 = dy2*dy2+dx2*dx2;
%第三通道
dy3 = im(y+1,x,3)-im(y-1,x,3);
dx3 = im(y,x+1,3)-im(y,x-1,3);
v3 = dy3*dy3+dx3*dx3;

%选择最强的梯度
if (v2 > v)
v = v2;
dx = dx2;
dy = dy2;
end
if (v3 > v)
v = v3;
dx = dx3;
dy = dy3;
end

%对应到18个方向中的一个
best_dot = 0;
best_o = 1;
for o = 1:9
dot = uu(o)*dx + vv(o)*dy;
if dot>best_dot
best_dot = dot;
best_o = o;
else
if -dot>best_dot
best_dot = -dot;
best_o = o+9;
end
end
end

%使用线性差值加到像素周围的四个直方图中
xp = (x-0.5)/sbin + 0.5;
yp = (y-0.5)/sbin + 0.5;
ixp = floor(xp);
iyp = floor(yp);
vx0 = xp-ixp;
vy0 = yp-iyp;
vx1 = 1.0-vx0;
vy1 = 1.0-vy0;
v = sqrt(double(v));
%左上角
if ixp >= 1 && iyp >= 1
hist(iyp,ixp,best_o) = hist(iyp,ixp,best_o)+vx1*vy1*v;
end
%右上角
if ixp < blocks(2) && iyp >= 1
hist(iyp,ixp+1,best_o)  = hist(iyp,ixp+1,best_o)+vx0*vy1*v;
end
%左下角
if ixp >= 1 && iyp < blocks(1)
hist(iyp+1,ixp,best_o) =  hist(iyp+1,ixp,best_o)+vx1*vy0*v;
end
%右下角
if ixp < blocks(2) && iyp < blocks(1)
hist(iyp+1,ixp+1,best_o) =  hist(iyp+1,ixp+1,best_o)+vx0*vy0*v;
end

end
end
%通过计算梯度方向上的和计算每个block中的能量
for o = 1:9
for x=1:blocks(2)
for y=1:blocks(1)
norm(y,x) = norm(y,x) + (hist(y,x,o)+hist(y,x,o+9))*(hist(y,x,o)+hist(y,x,o+9));
end
end
end

eps = 0.0001;
%计算特征
for x = 1:out(2)
for y = 1:out(1)
%右下角
n1 = 1/sqrt(norm(y+1,x+1)+norm(y+2,x+1)+norm(y+1,x+2)+norm(y+2,x+2)+eps);
%右边
n2 = 1/sqrt(norm(y,x+1)+norm(y+1,x+1)+norm(y,x+2)+norm(y+1,x+2)+eps);
%下边
n3 = 1/sqrt(norm(y+1,x)+norm(y+2,x)+norm(y+1,x+1)+norm(y+2,x+1)+eps);
%自己
n4 = 1/sqrt(norm(y,x)+norm(y+1,x)+norm(y,x+1)+norm(y+1,x+1)+eps);

t1 = 0;
t2 = 0;
t3 = 0;
t4 = 0;

%对比度敏感特征
for o = 1:18
h1 = min(hist(y+1,x+1,o)*n1,0.2);%截短
h2 = min(hist(y+1,x+1,o)*n2,0.2);
h3 = min(hist(y+1,x+1,o)*n3,0.2);
h4 = min(hist(y+1,x+1,o)*n4,0.2);
feat(y,x,o) = 0.5 * (h1+h2+h3+h4);%求和
t1 = t1+h1;
t2 = t2+h2;
t3 = t3+h3;
t4 = t4+h4;
end

%对比度不敏感特征
for o = 1:9
sum = hist(y+1,x+1,o) + hist(y+1,x+1,o+9);
h1 = min(sum*n1, 0.2);
h2 = min(sum*n2, 0.2);
h3 = min(sum*n3, 0.2);
h4 = min(sum*n4, 0.2);
feat(y,x,o+18) = 0.5 * (h1+h2+h3+h4);
end

%纹理特征
feat(y,x,28) = 0.2357*t1;
feat(y,x,29) = 0.2357*t2;
feat(y,x,30) = 0.2357*t3;
feat(y,x,31) = 0.2357*t4;
end
end

%%%%%%%%%下面开始视觉化的部分%%%%%%%%%%%%%%%
%对应于原程序的HOGpicture
bs = 20;
w = feat(:,:,1:9);
scale = max(max(w(:)),max(-w(:)));

% 给每一个方向做一个“条纹沟”
bim1 = zeros(bs, bs);
bim1(:,round(bs/2):round(bs/2)+1) = 1;
bim = zeros([size(bim1) 9]);
%bim(:,:,1) = bim1;
for i = 2:9,
bim(:,:,i) = imrotate(bim1, -(i-1)*20, 'crop');
end

% 通过将权重化的条纹加起来做成正权重的图案
s = size(w);
w(w < 0) = 0;
im = zeros(bs*s(1), bs*s(2));
for i = 1:s(1),
iis = (i-1)*bs+1:i*bs;
for j = 1:s(2),
jjs = (j-1)*bs+1:j*bs;
for k = 1:9,
im(iis,jjs) = im(iis,jjs) + bim(:,:,k) * w(i,j,k);
end
end
end
pos = im * 255/scale;

% %通过将权重化的条纹加起来做成负权重的图案
% w = -w;
% s = size(w);
% w(w < 0) = 0;
% im = zeros(bs*s(1), bs*s(2));
% for i = 1:s(1),
%   iis = (i-1)*bs+1:i*bs;
%   for j = 1:s(2),
%     jjs = (j-1)*bs+1:j*bs;
%     for k = 1:9,
%       im(iis,jjs) = im(iis,jjs) + bim(:,:,k) * w(i,j,k);
%     end
%   end
% end
% neg = im * 255/scale;

% 将图像加起来画一下
buff = 10;
pos = padarray(pos, [buff buff], 128, 'both');
%neg = padarray(neg, [buff buff], 128, 'both');
%im = uint8([pos; neg]);
im = uint8(pos);
clf;
imagesc(im);
colormap gray;
axis equal;
axis off;

toc

测试时所用的图像为



求得的HOG可视化以后为：



可以看到有很多竖线，说明竖直方向的HOG特征比较强，但总觉得有问题。将竖直方向的HOG特征分量置为零以后，即将可视化代码中的bim(:,:,1) = bim1;这句注销掉，

会得到下面的图像。



看着更合理一些吧。代码我是完全按照voc中的features.cc编写的，如果有谁发现有不一样的地方，请一定不吝赐教。下一篇应该会是latant-svm的matlab版本。