上面三篇高度相似图片检测没有讲到的一个方法就是基于图片特征点这个feature做图片相似检测。 图片特征点源于1999年Lowe 那篇”Object recognition from local scale-invariant features” 该文可以说是在图像识别形成了一前一后的格局，几年之后各个基于Lowe特征点的方法出来了，而原来基于颜色或者texture的方法就相对于naive了。以后的几年各种优化方法都出来了。 Lowe的方法是Scale-invariant feature transform 简称SIFT，从名字上来看这个方法对于图片的大小变化能起到很好忽略。这里有一点就是，这些特征点的主要目的就是 两张类似图片上的类似物体的特征点值是类似的，而不相似物体区域的特征点值是不一样的。在检测特征点优劣的时候主要就是检测一个repeatability。 Repeatability 就是只类似的物体特征点， 在经历不同的transformation之后的取回率，当然取回率高不一定好，还要看false positive，就类似信息检索。 transformation有很多种: 大小，旋转，颜色亮度，affine转变，occlusion什么的。 各种方法在针对不同的变换有不同的优化。

市面上有很多图片特征检测器和变种，下面是我用过的一些。

• Canny: 1986年出来的边检测器。现在还一直用，但不属于特征点范畴，更多用于基于boundary的方法。
• Harris corner detector: 边角检测器
• SIFT: 虽然这么叫，其实是通过高斯差获取scale space的最大最小值的位置。
• SURF: 类似SIFT，速度更快点。
• Harris-Affine 和 Hessian Affine: Mikolajczyk 04年的state of the art 方法。

检测器检测出来的特征点要去’形容’它，把他转变成特征向量。最长用的就是SIFT descriptor，不局限于sift检测器，上述方法都可以用它，只要知道 特征点的位置，大小。一般检测特征点，会在不同的scale space上，每个scale space出来的特征点大小会不一样。SIFT特征向量是一个维度128的向量，再加上一个 特征区域的角度。 如下图所示，中心点周围的有4个1616的格子，每个格子里有16个44的格子，会计算每个小格子的方向，划分成8个方向，每个4x4的小个子会成成一个长度为8 的向量，那总共有16个这样4*4大小的格子，所以最后会合成一个维度为128的向量。

通过获得sift向量。一张图片能很好的给压缩成几k的向量矩阵。在检测相似图片的时候，对两张图片的sift向量逐一比对，如果距离在特定阀值之下，就能确定是相似的图片。 但是这样的话计算量是很大的， 因为通常一张1200*960的图片会有700-2000的sift特征点。那如果我的数据库里有1000张图片要比对的话，那就是很大个开销， 而且，sift向量逐一比对的时候，需要计算每对之间的距离。所以很多时候，会通过bag of words 把sift向量分类成x个cluster (x小于sift向量数), 然后把一张图片通过 它在每个cluster所拥有的sift向量数目来形容。这样一张图片就变成了一个长度为x的直方图。这样比对起来就会快了。虽然这样也会有很多信息丢失。但在实践中，这招非常有用。

基于特征点的方法的有点是对于旋转，affine上的变换会比前几个方法容忍度更高，在尺寸，亮度上的容忍也非常棒, 而且如果出现物体的遮盖occlusion，有时候也能检测出来。就是可能速度上没有前者的快。 这个特征点的方法不错，不过有一个致命的缺点就是所谓的semantic gap，往往你有很多你有很多sift的cluster，但是他们的空间信息都是丢失掉的。并且，一个sift和另一个sift之间的相对 距离也不会给利用到。会出现false positive，就是图片和另外一个外形上外向不相似的图片的相似度很高。

家里没网，只好扯扯淡了。

虽然一个是无监督学习，一个是监督学习，但是我感觉两个算法是近亲。 knn的就是找最近友邻，通过数据点和数据点之间的位置关系，而随机森林的话，是通过信息熵把类似的数据都放在一个叶子上。 其实这个感觉以前看到在集智俱乐部(其实就去过一次)有一面之缘的大牛的豆瓣博客上写到的，但是当时没有领悟到里面的缘由。毕业论文用到的分类算法就这两个，现在就觉着他们是近亲。

两个算法各有优缺点，先说knn吧，不需要训练，只要有个好点的数据结构储存数据让查询的效率高一些就好(以前的博客有写)，找到的总是相似度最高的友邻。这样的好处就是不”吃亏”啦，没有loss。唯一要确定的的参数也就是取多少友邻介一个，这种naive算法给像我这种初学机器学习吐阳图性破的人用再好不过了。不过在测试的时候的速度是knn的软肋。

那随机森林呢，你要随个鸡鸡森林的话，设置起来就难啦，首先森林是基于树木的，你要学会ID3, C5(不是炸弹)决策树什么的，信息熵，gini啊什么的,各种数据分割的criteria，麻烦的很，要设置参数也多，什么多少层熟，一个森林多少树等等。不过参数多也有好处，就是可适性高。随机森林的话，有点在于速度超级快，从root到leaf的这段过程就是一个knn最近友邻的查找过程。一般在每个节点就是一个很简单的比较，所以速度超级快。结果会有loss，不一定是最相似的友邻。

介绍完两个算法的脾气后，得说说我最近几天和他们相处的感觉。knn适用于数据少的时候。RF用于数据多的时候。首先速度上来说knn就是这样的，另外就是，数据少的时候，友邻质量对你来说特别重要，你都没几个朋友了还都是帮土匪的话那你也就完了。RF作为一种boosting方法，本身就是个statistic概念很强的东西，通过对每个树加随机因子(数据上的或者特征上的随机，也有别的随机方式，譬如增加噪音什么的)让每棵看着差不多，用起来不一样:decorrelation。 这样通过增加variance的方法让结果更好点。但是如果数据少的话，就会适得其反啦。类似的，数据多的时候，knn拿到的优质友邻都很类似，那其实结果就是拿到了更少的友邻，换个角度看的话，你在射交网络天天关注的就几个人，这几个人还都是吃喝拉撒在一个地方的近亲，那他们给的信息大部分都差不多的，这时候你需要weak tie的友邻啦(好吧，我又扯到射交网络了，sorry)。而这个knn的问题rf就能弥补。这里你要说了，你不是说友邻质量很重要么？对，那是在数据少的时候，因为少，从概率上讲友邻的他们之间相似程度不会很高。

蛋疼无极限，想在c++ 里跑octave代码，找到了两个帖子，拼拼凑凑终于可以跑了！ 帖子1 还是邮件列表里的私藏。 照着上面的代码跑一下会发现没有装octave的头文件

    Hi,

    I tried to find information on how to call an octave .m file from a  
    C++ program. There is some information around, but most of it is no  
    longer up-to-date, it was not trivial to get it working. In this post  
    I want to explain how to do it, maybe it is useful to somebody.

    I have a file called exampleOctaveFunction.m which looks like this:

    -------------------------------------
    function [resultScalar, resultString, resultMatrix] =  
    exampleOctaveFunction (inScalar, inString, inMatrix)

    resultScalar = (inScalar * pi);
    resultString = strcat ('Good morning Mr. ', inString);
    resultMatrix = (inMatrix + 1);

    endfunction
    -------------------------------------

    I have a file called how-to-call-octave.cpp which is this:

    -------------------------------------
    #include <octave/oct.h>
    #include <octave/octave.h>
    #include <octave/parse.h>
    #include <octave/toplev.h> /* do_octave_atexit */

    int main (const int argc, char ** argv)
    {
    const char * argvv [] = {"" /* name of program, not relevant */,  "--silent"};

        octave_main (2, (char **) argvv, true /* embedded */);

        octave_value_list functionArguments;

        functionArguments (0) = 2;
        functionArguments (1) = "D. Humble";

        Matrix inMatrix (2, 3);

        inMatrix (0, 0) = 10;
        inMatrix (0, 1) = 9;
        inMatrix (0, 2) = 8;
        inMatrix (1, 0) = 7;
        inMatrix (1, 1) = 6;
        functionArguments (2) = inMatrix;

        const octave_value_list result = feval ("exampleOctaveFunction",  
        functionArguments, 1);

        std::cout << "resultScalar is " << result (0).scalar_value () << std::endl;
        std::cout << "resultString is " << result (1).string_value () << std::endl;
        std::cout << "resultMatrix is\n" << result (2).matrix_value ();

        do_octave_atexit ();
        }
    -------------------------------------

    And a little readme file, called readme.sh, that explains how to  
    compile and run this simple example:

    -------------------------------------
    #! /bin/bash

    #
    # Make sure you have octave installed. On my system, Ubuntu 8.10, I installed
    # the packages octave3.0 and octave3.0-headers.
    #
    # To compile, type
    #
    make how-to-call-octave.o
    #
    # which, on my system, is equivalent to
    #
    #   g++ -c -o how-to-call-octave.o how-to-call-octave.cpp
    #
    # To link, type
    #
    g++ -L /usr/lib/octave-3.0.1/ -l octinterp -o how-to-call-octave  
    how-to-call-octave.o
    #
    # Adjust the directory in the -L directive according to the configuration of
    # your machine.
    #
    # To run, type
    #
    ./how-to-call-octave
    #
    # It should output something like
    #
    #   resultScalar is 6.28319
    #   resultString is Good morning Mr. D. Humble
    #   resultMatrix is
    #    11 10 9
    #    8 7 1
    #
    # Hope this is of help to somebody.
    #
    -------------------------------------

    You can also execute the readme.sh file and it will compile and run for you.
    Good luck, Arjan.

apt-get 之后发现还是跑不起来, 在帖子2里找到了解决方法

I started out trying to run the octave example as a c++ file(i.e. no Qt) and received the same errors. The problem is you have to tell the linker where the libs are located. As root, run "ldconfig /path/to/octave/libs". This is not permanent because the next time ldconfig is run it will not find the octave libs. 

To make sure the cache gets updated every time ldconfig is run (Ubuntu is a Debian derivative so I'm assuming the procedure is the same) do this:


sudo -i
cd /etc/ld.so.conf.d
touch octave.conf
vi (your fav. editor) octave.conf 
add a line: /path/to/octavelibs
save file
run ldconfig

设置后就能跑啦。最近的博客越来越水了。我又要烂掉了。

图片中的物体探测一直是比较热的课题，通过物体的检测可以对图片自动生成标签以及分类，可以很好的用在图片检索的领域。在做的毕业设计(也是论文)就涉及到这方面的。今天看一篇叫Class-Specific Hough Forests for Object Detection的state of the art，发现原来导师叫我用到的方法和这个差不多，只是增加了某些小trick用于对特定的数据集优化。

Hough Forests。 其实就是一个random forests 也就是随机森林,每条数据是一小块图片的description(patch) 树的每个节点对数据的分割通过information gain，努力把相同类的patch和相似的patch分到同一个枝头上去。每次通过数据集,随机采集部分features用于分割，每个这些features的分割点生成一个pool。然后对每个feature的分割点做一个信息量的检测，就像一个标准的random forests一样， 这里信息量的计算通过一下程式(这么说比较台湾腔 夯): $$ Entropy({c_i}) = -c\cdot \log c -(1-c)\cdot \log(1-c) $$ 这里c就是patch是属于物体c的概率，1-c就是negative cases。o 简单吧，其实在原来的论文里还有一个东西，就是entropy只是衡量split的优差的其中一个feature，还有一个feature是offset也就是patch到中心点的距离，通过增加这个衡量点，相似距离的patch会在一起哦。测试的结果是，加上这个feature，precision会有很好的提高。 怎么计算这个feature，其实很简单，用Root mean square(RMS)把split后的patch到中心点的距离计算下，split之后小的RMS更优。到达叶子的数据，每个patch会有一次投票的机会给出中心点的位置。

训练之后，每个测试图片会被分成不同的patch，通过hough forests 到达叶子的patch生成不同中心点的位置(根据在那叶子里的中心点patch). 最后，通过投票，某个物体的中心点的票数会很高，就这样我们就能找到一个物体的中心点了。 如果一个patch是不是物体的一部分，生成的中心点是错误的，但是中心点的位置其实是随机的，也就是说如果有很多很多非物体patch生成的map是没有一个明显的中心点的。而如果有很多物体中某部位的patch，生成出来的中心点是有一个很大的偏向性的，往往会往某一个中心聚拢。

这里有一个问题就是怎么计算一个patch到中心点的位置。首先训练数据一定得是每张图片只有一个物体。每个patch计算自己中心到物体的中心并保存x轴和y轴的移动。

道理很简单，给我的代码。 opencv这个库让我又爱又恨啊，很好用，很多很复杂的算法只要一行代码就ok了，但是python binding里有很多bug。而且很多错误你不一行行代码拆下来是不会看到的。。。 首先，不想使用论文里的hough forest，因为感觉上在opencv里你得自己hack下，不值得。第二貌似没找到patch的descriptor，所以在自己的实验里，只用到了特征点的descriptors，像sift和surf。使用这个的有点就是，可以做到orientation insensitive和scale insensitive。 其实，论文里的方法是没有能做到scale insensitive的，因为距离是个relative的东西，论文里，通过使用不同的scale进行测试，获得最优的结果作为最终结果。 而view insensitive更没有解决方法，按他们的方法，只有通过增加不同view的训练数据来弥补吧。

但是！导师教给我的方法能避免这个问题。首先sift和surf本身是view insensitive的，而在计算距离的时候通过对特征点本身的大小对距离做个归一化，这样就能让系统获得对不同大小的物体的识别能力。 另外一个trick就是sift之类的特征点自己会带有gradient orientation,也就是梯度方向，这个信息是在计算距离的时候必须使用到的。因为相同descriptor的特征点可能会有不同的梯度方向，如果你只使用了图片原始的x和y轴保存距离的话，预测的中心点会偏离本来的实际位置。下面的图片就是个case:

我有一个圆圈，底部是白色的，如果是sift当特征点的检测的话，每个圆圈边上的点的特征值是一样的，而他们的梯度方向可能会不一样，如果不使用梯度方向的话，a作为训练的数据，而比作为测试数据，中心点预测的结果会在圆圈的外面(细线部分)。 如果要避免以上问题，只有通过rotation matrix 让每个特征点更具自己的梯度方向生成自己的x和y轴。 并且每个特征点的相对于中心点的距离要在新的x和y轴上进行计算。

• $$x_{new} = x \cdot cos(angle) - y \cdot sin(angle)$$
• $$y_{new} = x \cdot sin(angle) + y \cdot cos(angle)$$

如果需要从新轴返回到原始的x y轴只需要把angle乘以-1就可以了, 上面的angle是以radian为单位的。

假设上图是更具特征点的梯度方向做的x和y轴的转换，转换之后，可以发现点离圆圈的中心的距离不变，但是在x轴和y轴的移动发生了变化。

下面是测试方法的代码，包括训练和测试。数据集的话是使用ucic里的汽车数据，训练集里有positive cases 和negative cases，这里我只用到了positive cases。也没有是用到hough forest做训练，只用到了最近友邻。 每次选择最相近的友邻做出投票。在训练数据里，所有的车子中心点都在图的当中，因为每张图片里的车子几乎撑满整张图，而所有图片的尺寸都是40x100，在训练数据里，有两部分，一部分是没有scale变化的车子数据库，而另外一部分里车子的大小会有变化，用于测试方法对于scale变化的鲁棒性。

下面是python代码

import cv2
import math
import sys
import os
import numpy as np
kernel=np.array([[0,0,0,0,0,0,0,0.0002,0.0009,0.0014,0.0016,0.0014,0.0009,0.0002,0,0,0,0,0,0,0],[0,0,0,0,0,0.0005,0.0021,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0021,0.0005,0,0,0,0,0],[0,0,0,0.0000,0.0016,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0016,0.0000,0,0,0],[0,0,0.0000,0.0020,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0020,0.0000,0,0],[0,0,0.0016,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0016,0,0],[0,0.0005,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0005,0],[0,0.0021,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0021,0],[0.0002,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0002],[0.0009,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0009],[0.0014,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0014,],[0.0016,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0016],[0.0014,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0014],[0.0009,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0009],[0.0002,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0002],[0,0.0021,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0021,0],[0,0.0005,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0005,0,],[0,0,0.0016,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0016,0,0],[0,0,0.0000,0.0020,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0020,0.0000,0,0],[0,0,0,0.0000,0.0016,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0016,0.0000,0,0,0],[0,0,0,0,0,0.0005,0.0021,0.0031,0.0032,0.0032,0.0032,0.0032,0.0032,0.0031,0.0021,0.0005,0,0,0,0,0],[0,0,0,0,0,0,0,0.0002,0.0009,0.0014,0.0016,0.0014,0.0009,0.0002,0,0,0,0,0,0,0]])

FLANN_INDEX_KDTREE = 1  # bug: flann enums are missing
AUTOTUNED = 255
flann_params = dict(algorithm = FLANN_INDEX_KDTREE, trees = 4)


def get_key_points_from_img(img):
    if isinstance(img,str):
        img = cv2.imread(img)
    im = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

    surf = cv2.SURF(1000, 4, 2, True)
    keypoints, descriptors = surf.detect(im, None, False)
    l_blue = cv2.cv.RGB(200, 200, 250)
    descriptors.shape = (-1, surf.descriptorSize())

    # y x
    center = (im.shape[0]/2,im.shape[1]/2)
    attributes = []

    for i in range(len(keypoints)):
        kp = keypoints[i]
        # (x y)
        scale = kp.size
        #cv2.circle(im, (int(kp.pt[0]), int(kp.pt[1])), int(kp.size/2), l_blue, 2, cv2.CV_AA)
        gradient = math.radians(kp.angle)

        #calculate the distances x and y in the new axis(rotated by the gradient)
        d_x = center[1]-kp.pt[0]
        d_y = center[0]-kp.pt[1]
        tx = math.cos(gradient) * d_x - math.sin(gradient) * d_y
        tx /= scale
        ty = math.sin(gradient) * d_x + math.cos(gradient) * d_y
        ty /= scale
        # generate data for further experiments
        attributes.append((kp.pt[0], kp.pt[1], scale, gradient, ty, tx, 0, 0, 0, 0, 1))
    return (descriptors,attributes)



if __name__ == "__main__":

    dir_path = "./CarData/TrainImages/"
    files = os.listdir(dir_path)

    #build train dataset
    train_descriptors = []
    train_attributes = []
    for f in files:
        if f.startswith("pos"):
            (descriptors, attributes) = get_key_points_from_img(dir_path+f)  
            train_descriptors.extend(descriptors) 
            train_attributes.extend(attributes)

    print len(train_descriptors)
    train_descriptors = np.array(train_descriptors)
    #build pnn model
    flann_index = cv2.flann_Index(train_descriptors, flann_params)

    dir_path = "./CarData/TrainImages/"
    dir_path = "./CarData/TestImages/"
    dir_path = "./CarData/TestImages_Scale/"

    files = os.listdir(dir_path)
    count1 = 0
    count2 = 0
    for f in files:
        img = cv2.imread(dir_path+f)
        (descriptors, attributes) = get_key_points_from_img(img)  
        im = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
        mapa = np.zeros(im.shape, dtype=float)
        indexes, dist = flann_index.knnSearch(descriptors, 15, params={})  
        for i in range(len(descriptors)):
            for j in range(len(indexes[i])):
                scale = attributes[i][2]
                gradient = attributes[i][3]
                d_y = train_attributes[indexes[i][j]][4]*scale
                d_x = train_attributes[indexes[i][j]][5]*scale
                tx = math.cos(-gradient) * d_x - math.sin(-gradient) * d_y
                ty = math.sin(-gradient) * d_x + math.cos(-gradient) * d_y
                predicted_center_x = int(round(attributes[i][0]+tx))
                predicted_center_y = int(round(attributes[i][1]+ty))
                if (predicted_center_x>=0) and (predicted_center_y>=0) and (predicted_center_x<mapa.shape[1]) and (predicted_center_y<mapa.shape[0]):
                    mapa[predicted_center_y][predicted_center_x] += 1
                    count1 += 1
                else:
                    count2 += 1

        #apply a matlab disk filter to the output map for visualization
        mapa2 = cv2.filter2D(mapa,-1,kernel)*5
        cv2.imshow("original",im)
        cv2.imshow("image", mapa2)
        cv2.waitKey()

一大坨kernel变量，其实就是matlab的disk filter: fspecial(‘disk’, 20) opencv里就没有合适的filter可以用来平滑图片的!。 这里使用到了flann用来做最近友邻的模型，由于数据不是特别大，所以这里也没什么优化。在训练的时候，会把每个特征点距离到中心点的距离计算一下，计算完会根据特征点的大小把距离缩放一下，而当在测试的时候也会根据每个特征点的大小把距离放大。这样就能达到scale insensitive的目的了。 每个特征点会根据友邻做投票，投票的结果会映射到一个和图片一样大小的map上面。 通过filter2d过滤之后，会出现一大坨模糊的东西，哪里最亮哪里的投票数就最多，是中心点的可能性也最大。

下面就是一些实验图片的结果。

貌似结果很好的，但是其实不是的，因为这里没用到false cases，会出现很多false positive，而且，由于是使用sift detector, 本身控制不了哪些特征点是需要的哪些是不需要的，在测试后的时候，如果测试图片的某个非物体部位的特征点很多的话(细节比较多)，会产生很多假票, 像最后一张图的两根柱子。在下次的试验力，一定得吧那些false cases包含到训练库里，当测试数据碰到的一个false case的票，直接把权值归0，这样的话应该能对结果有很大的改观。还有就是这里，每张投票的权值都是1，可以根据友邻的距离设置投票的权值，近的投票全职大，这样结果应该会好一些。

opencv用起来真的没有matlab舒服，就连一个类似matlab的imregionalmax函数都没有，这个函数可以计算local maximum，这样的话，选取峰值最高的几个点作为预测结果会有很理想的结果。所以这套代码在这里也不能用来获得和真实中心点做比较。

最近需要使用最近友邻对高纬度的数据进行查询。在实验的过程中发现相对于高纬度的数据，kd树的查找速度其实没有线性查找来的快的，但是在大数据面前，线性查找的开销是很大的。我在网上找了一些相关资料，最后发现有个叫flann的库很不错，这个库可以让根据用户设置的查找准确性确定性和用户提供的数据自动生成一个合理高效的数据结构，让每次搜索的速度比线性查找快几个数量级，同时丢失一些准确度。程序上自动的对参数和数据结构的设置，让使用者不用担心该使用什么样的方法来解决问题。

网上有各类优化的最近友邻查找方法的论文，很多时候针对每个问题，用的方法不一样能产生不同的优化。通过论文实验报告，有2个方法是在比较不错的:

• Randomized kd-tree algorithm: 不同于普通的kd-tree，这个方法是通过生成一组(大于等于1)的kd树来储存数据。该方法是通过在最之前的D次随机数据分割从而获得方差最大的分割(D=5)。当搜索的时候所有树会根据和数据的距离进行排列，距离近的树先被排查，只会查找固定数目的leaf。
• Hierarchical k-means tree algorithm: 这个算法就是每次把数据进行k组分割(k-means) 每分割完的小组继续进行k组分割，直至每组数据小于K。这里随机部分是:首先查询数据会在树里过一遍，到达叶子的时候会生成一个优先队列，在队列里，查询数据到达最有叶子里的路径里的所有节点里没有访问过的branch, 队列的排列方式根据他们的中心位置和查询数据的距离从小到达排列。每次会从优先队列里提取一个节点从那里一直到树叶进行查找最近友邻，并同时在优先队列里插入新的未访问branch。查找的次数根据用户的提供的准确率百分比。

最有参数的自动选择没有怎么整明白(刻意关系)，貌似是使用10%的数据做一个检测，参数的选择会在固定的几个之中进行选择。用户使用的时候往往需要提供一些使用场景的参数，比如: 查找速度，建表速度，内存使用速度，更具这几个参数，他会自动给你找到最优的方法和方法参数。只提供使用场景的这些参数更直观的让不是怎么了解算法本身的人能够用最合适的算法做最合适的事了。具体cost计算方程在论文的程式1中。

使用的感脚良好，很简单，速度也靠谱。有python matlab 等的binding，opencv里还自带了这个库！在mac上安装matlab上安装的时候碰到几个小问题。安装步骤可以参考http://www.cs.uky.edu/~jacobs/tips/flann_matlab.html。当中需要装一个patch让matlab里的mex可以能编译。地址在这里装完之后更新下配置文件就能在matlab上编译c++文件了！

这个库的使用也是很容易的:

[index, search_params, speedup] = flann_build_index(mat_train_data(:,:), struct('algorithm','autotuned', 'target_precision',0.95, 'build_weight',0.01, 'memory_weight',0))
%这里表示无所谓内存开销，建立索引开销重视一般。
%搜索        
[result, ndists] = flann_search(index, mat_train_data(:,1), 5, search_params);
%保存
flann_save_index(index,'train.idx')

不过还是有很蛋疼的地方，比如，有时候你得save一下再重新load之后才能搜索, 还有就是里面数据和matlab的数据表现形式不一样啊，让我废了好久才弄明白，一般数据分析总是会把数据一行行保存，每个row是一条数据。但是这里他是反人类的反过来了！使用的时候要非常注意。

参考:

• 1 FAST APPROXIMATE NEAREST NEIGHBORS WITH AUTOMATIC ALGORITHM CONFIGURATION
• 2 http://www.cs.uky.edu/~jacobs/tips/flann_matlab.html
• 3 http://www.cs.ubc.ca/~mariusm/index.php/FLANN/FLANN

以前就想做一个马赛克合成的图片了，最近不务正业，就写了一个 方法很简单:

• 抓取一组图片的并缩放到统一大小
• 计算每张图片的rgb值的中值
• 所有的图片中值生成kd-tree
• 加载目标图片，并缩放到想要的结果图片的大小，然后创建结果图片(np.zeros)
• 每一小格结果图片的大小和抓取的小图大小一样,计算目标图片每小格的rgb 中值
• 查找最想进的小图，并用小图的rgb值填充到结果图片里
• 生成结果图片

马赛克图片生成代码:

import Image
import os
import numpy as np
from scipy import spatial

#calculate a mean rgb of an image
def meanRGB(img):
    rgb = np.array(img.getdata())
        mean_r = np.mean(rgb[:,0])
        mean_g = np.mean(rgb[:,1]) 
        mean_b = np.mean(rgb[:,2])
        return [mean_r, mean_g, mean_b]

if __name__ == "__main__":
    #grid size
    mosaic_img_size = (20, 20)
    #size of target image to generated
    target_img_size = (1000, 1000)
    target_img = Image.open("./VanGogh-starry_night_ballance1.jpg")
    target_img = target_img.resize(target_img_size)
    #grid images directory
    directory = "./icons/"
    listing = os.listdir(directory)
    img_list = []
    color_list = []
    #load images and calculate mean rgb
    for infile in listing:
        if infile.endswith("jpg"):
            img = Image.open(directory+infile)
            img = img.resize(mosaic_img_size, Image.ANTIALIAS)
            img_list.append(img)
            color_list.append(meanRGB(img))

    #create kdt for future pnn query
    data = np.array(color_list)
    tree = spatial.KDTree(data, leafsize=5)

    count_x = 0
    count_y = 0
    new_image = np.zeros((target_img_size[0],target_img_size[1],3), dtype=np.uint8)
    target_img_data = np.reshape(np.array(target_img.getdata()), (target_img_size[0],target_img_size[1],3))
    length = mosaic_img_size[0]*mosaic_img_size[1]

    #calculate mean rgb of each grid in the target image, search a nearest neighbor and fill the new image with
    #the rgb values of the grid image
    while count_x+mosaic_img_size[0] < target_img_size[0]:
        count_y = 0
        while count_y+mosaic_img_size[1] < target_img_size[1]:
            mean_r = np.mean(target_img_data[count_x:count_x+mosaic_img_size[0],count_y:count_y+mosaic_img_size[1],0])
            mean_g = np.mean(target_img_data[count_x:count_x+mosaic_img_size[0],count_y:count_y+mosaic_img_size[1],1])
            mean_b = np.mean(target_img_data[count_x:count_x+mosaic_img_size[0],count_y:count_y+mosaic_img_size[1],2])
            img_index = tree.query(np.array([mean_r,mean_g,mean_b]))[1]
            new_image[count_x:count_x+mosaic_img_size[0],count_y:count_y+mosaic_img_size[1],:] = np.reshape(np.array(img_list[img_index].getdata()), (mosaic_img_size[0],mosaic_img_size[1],3))

            count_y += mosaic_img_size[1]
        count_x += mosaic_img_size[0]
    
    #image created
    im = Image.fromarray(new_image)
    im.show()
    im.save("./1VanGogh-starry_night_ballance1.jpg")

里面有些小bug，但不影响使用，譬如最右最下有一圈黑色的，是代码中某个逻辑没写好，没有填充进去，但是不管它了- -(貌似把小于号改成小于等于就可以了)

原图: 效果: 768的向量空间长度，图片量一定要大一些，要不然的结果会是，颜色很单一，拼出来的马赛克效果不好。自己的友邻太少了，可能结果会不理想，所以这里我抓了豆瓣活动上的豆友，570多位。。

最后一个要介绍的图片指纹生成方法: pHash, 这里是pHash算法简要:

• 图片指纹灰值化
• 图片缩放到32x32
• 计算机视觉每个图片的Discrete Cosine Transform (DCT, type II)。
• 只保留最左上部的8x8的DCT系数(有32x32大)，这样保存了频率最低的那部分(图片信息的大部分)
• 计算中值
• 生成64维2维码，0表示系数小于中值，反之为1.

说到DCT，其实就想傅里叶转换一样的，把信号转换成很多不同频率和振幅的正玄曲线相加的结果。但是DCT只是用cosine函数。这样的话图片高频率部分会给擦去，只保留低频率部分。因为在给图片操作的时候，往往低频率的DCT系数会给保留，而且大部分图片的信息会给保留在这些低频率DCT系数里。(JPEG压缩也用这个方法来对图片进行压缩)。 DCT 会生成8x8的系数表，那最左上方的表示最低频率的元素，也是最重要的，越往右下的表示相对频率稍高的元素(好多信号出来的东西啊，头好大) 反正到最后能生成一个图片指纹，从别的作者的实验结果来看，效果是相对较好的。在这个网站上phash.org能下到开源的代码，在网站的demo上可以可以试试看计算图片的相似性，通过测试，只有DCT比较靠谱，另外两个的结果非常糟糕。相同的，这个算法的信息量为2^64，冲撞率应该不高，吧？

phash.org网站上是用c/c++写的代码，没有python的binding，在https://github.com/polachok/py-phash/可以找到一个简陋的python binding，里面有我们需要的DCT的计算方法，调用一个函数就能直接获得指纹数值。测试了几下，效果还不错。不过还是需要很大量的图片用语测试false positive的概率。 在github上能看到作者写的使用方法说明，我这里具体的代码如下:

import pHash
import sys

if __name__ == "__main__":
    hash1 = pHash.imagehash(sys.argv[1])
    hash2 = pHash.imagehash(sys.argv[2])
    print 'Hamming distance: %d (%08x / %08x)' % ( pHash.hamming_distance( hash1, hash2 ), hash1, hash2 )

图片缩小之后其实保留个部分图片的信息，可以用来当做图片的指纹，在查找相关资料的时候发现了一个很不错的博客 在这个博客里，几乎所有的和图片相似性计算的算法都在了(除了基于语义的)。超级不错。里面的很多算法的实现是用pearl的，而本身抱着学习的态度我决定重新用python对其中的某些算法实现一遍。

这里是图片缩小指纹的算法:

• 缩放到160x160的图片
• 灰值化
• 模糊图片，把噪音去除掉
• equalize，让图片反差更大
• 计算图片平均值
• 缩放到16x16的图，并且降到2维，更具图片pixel值:大于平均值的为1，反之为0

这样就生成了图片的指纹，信息总量有2^256, 哈希冲撞应该不大。 计算的图片相似度的时候直接使用xor。公式可以是这样: 1 - xor(p1,p2)/256

具体python实现代码如下:

from PIL import Image
import sys
import ImageFilter
import ImageOps 
import numpy as np

def normalize(img, size = (160, 160)):
    return ImageOps.equalize(img.resize(size).convert('L').filter(ImageFilter.BLUR))

def calculate_fingerprint(img):
    img = normalize(img)
    mean_value = np.mean(np.array(img.getdata()))
    img = img.resize((16,16))
    return np.array(img.getdata())>mean_value

def calculate_simility_by_path(img_path1, img_path2):
    img1 = Image.open(img_path1)
    img2 = Image.open(img_path2)
    fingerprint1 = calculate_fingerprint(img1) 
    fingerprint2 = calculate_fingerprint(img2) 
    return 1 - (0.0+np.sum(np.logical_xor(fingerprint1, fingerprint2)))/fingerprint1.size

if __name__ == "__main__":
    print calculate_simility_by_path(sys.argv[1],sys.argv[2])

更好的implementation可以在这个链接找到，虽然做法上不一样，但是原理应该都差不多。在代码的注释里，能看到这个类算法的优点和缺点，如果图片相似度较高的话，识别起来没什么问题。最大的缺点是，很多false positive，如果图片是一张大海的话 会出现一下情况

1111111111111111
1111111111111111
1111111111111111
1111111111111111
1111111111111111
1111111111111111
1111111111111111
1111111111111111
0000000000000000
0000000000000000
0000000000000000
0000000000000000
0000000000000000
0000000000000000
0000000000000000
0000000000000000

作者还说到对于大规模数据的话，性能会降低很多。但是貌似这个能用bloomfilter解决，现在先不管他

对于颜色直方图差的比对能很好的检测出图片之间的色差变化，那对于两张环境，光线，角度长不错的图片，能很好的检测出他们的相似度。 这是一个很老的方法，简单粗暴，非常使用。原理在这篇博客里很清楚的讲解了。 算法很简单:

• 把两张图片缩放到统一的大小，这里可以是256x256， 并且统一成rgb模式。
• 统计每张图片的每个颜色的出现次数。那这里因为是rgb模式，有3*256=768 种颜色，生成直方图
• 对两张图片的直方图进行比对。如果统计结果一致的话，相似度+1

这里有个问题就是对空间信息的丢失，原作者通过把图片分割成16个小方格来解决这个问题。(split_image函数)

具体代码如下:

from PIL import Image
import sys
def normalize(img, size = (256, 256)):
    return img.resize(size).convert('RGB')

def sim_hist(hist1, hist2):
    assert len(hist1) == len(hist2)
    return sum(1 - (0 if hist1 == hist2 else float(abs(hist1 - hist2))/max(hist1, hist2)) for hist1, hist2 in zip(hist1, hist2))/len(hist1)


def calc_similar_by_path(path_img1, path_img2):
    li = normalize(Image.open(path_img1))
    ri = normalize(Image.open(path_img2))
    return sum(sim_hist(l.histogram(), r.histogram()) for l, r in zip(split_image(li), split_image(ri))) / 16.0



def split_image(img, part_size = (64, 64)):
    w, h = img.size
    pw, ph = part_size
    assert w % pw == h % ph == 0
    return [img.crop((i, j, i+pw, j+ph)).copy() \
        for i in xrange(0, w, pw) \
        for j in xrange(0, h, ph)]

if __name__ == "__main__":
    print calc_similar_by_path(sys.argv[1],sys.argv[2])

这里的fingerprint就是每个图片的直方图。 这个指纹的信息量的话应该挺大的(至少有2^768) 理论上来讲应该false positive的可能性很小，但是这个结论应该是错的，因为在图片里很多颜色是出现的概率是很小的，大部分颜色都集中于某些值中。

代码和原作者的差不多，实验结果是不错的，很简单，但是这个算法有个致命的问题就是对颜色的过分依赖。对于图片角度的变化的容忍度挺高，但是如果颜色或者图片光线变化稍大就不能很好的检测出相似度了，会把相似的图片判断为不同的。最简单的方法就是弄个灰度图，这方法基于rgb的，就无解了。还有个小问题就是，false positive，我觉得如果两张白底黑字的文字图片，不管字是多么不一样，但是他们的相似度应该很高(在有的应用上，这个结果是很希望的)

在db实习了2个多月很开心。老实说，在那里做的项目没有想象中的理想，觉得惭愧啊。前段时间清风老师有和我提起一个敏感图片检测的需求。当时候有试着用SIFT和min-hash实现，但是结果不是很理想。在公司和在家自己弄有很大的区别，时间上的紧迫，很多时候虽然头不会催，但是你看着别人在交东西，自己却没有，压力很大，很多时候就需要move on。结果就是这个项目没有完成, 感脚好对不住清风老师啊！真的不会再爱了！

这几天闲下来，我突然有想好好研究研究这方面的东西了，看了些paper，觉得有很多实现方法，结果有好的有坏的，我准备都试着去实现下，把他们的有点和缺点都记录下来，希望能找到一个很好的解决方案，如果运气好的话。

Near-duplicate image(高相似图片)的检测不比完全一样的图片检测来的简单，后者可以直接用哈希生成像MD5类似的指纹，然后保存每个图片的时候也保存那个指纹，这样在查找的时候只要比对指纹就可以了，这样的话速度上会有很大的提升。

那用什么样的办法能有效的比对图片的相似度呢? 方法有很多。 首先要说的一点是，在比对的时候，速度是非常重要的，所以，一般都是通过指纹(fingerprint)技术把一张图片合理的压缩成一个容量占用很小的方便计算相似度的数据集。 前面说过md5是不能用在这里的，为什么呢？因为一个微小的变化会是两个图片之间的MD5完全不一样。而在这里要做的是:

• 相同图片的指纹要一样
• 类似图片的指纹也要类似
• 完全不相同的图片指纹的差别很大

做到以上几点，那么图片的相似度的识别就完成了，但是要找到一个函数f做到以上这种方法很难。这也是这里要慢慢探索的。