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

来源：互联网发布：数据分析师职业收入编辑：程序博客网时间：2024/05/16 05:20

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

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

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

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

function feat = HOG(input,sbin)ticuu = [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                endend%通过计算梯度方向上的和计算每个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    endendeps = 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;    endend%%%%%%%%%下面开始视觉化的部分%%%%%%%%%%%%%%%%对应于原程序的HOGpicturebs = 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  endendpos = 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版本。

0 0