Posted on | In

问题

最近在安装python软件包的时候，经常会遇类似这样一个问题。比如对于ipython,机子本身安装的版本是1.2.1，显然太低，不足以跑jupyter，尝试着用pip安装，却发现下载的过程一路畅通，但是安装的时候却总是会报这样一个错误：

1
......'Not uninstalling ipython at /usr/lib/python2.7/dist-packages, owned by OS'......

最终显示了ipython已经成功安装，但是查看version的时候却仍然是1.2.1，弄得我十分头大。

原因

后来研究了一下，发现其实造成这个问题的原因很简单，就是因为ubuntu的apt-get版本与pip安装的版本冲突。由于apt-get的’地位‘要比pip高那么一点（毕竟亲爹），系统会优先使用apt-get 里面的软件包。但是apt-get 里的软件通常特别的老旧，完全无法跟得上python包的更新速度，pip虽然版本新，但是却不能删除apt-get 的老旧版本。。。。。。这就造成了错误中提到的的无法卸载的问题了。

解决

至于此，解决问题的方法已经很清楚了，只要手动卸载掉apt-get里对应的过时的软件包即可。

Posted on | In

jupyter其实就是ipython notebook的另一个版本，是一个很强大的基于ipython的python代码编辑器，python文档查看器。他可以部署在网页上，可以非常方便的对文件进行查看、下载，并且对python文件进行在线编译，甚至是远程连接。。。。。。用他编写的python文件本身就是一个强大的开发文档。更重要的是很多基于python开发的开源程序都热衷于用.ipynb格式的文件作为文档（比如caffe）。因此掌握jupyter的使用也尤为重要。

Ipython notebook

简介

首先我们来探讨下jupyter的前身，ipython notebook。直接输入命令$ipython notebook 即可在浏览器中自动打开一个界面，显示了当前路径下的文件树。在这当中我们可以自由的编辑、处理文件，编写运行python文件等等。我们一方面可以用它来编程，另一方面也可以把他当成一个远程文件管理器，界面好看，用起来也十分的方便。

既然如此，我们就试着在服务器上搭建一个ipython notebook服务。

服务器配置

详细的文档可见ipython官方文档，这里是结合文档进行的配置。

配置文件

配置文件默认在~/.ipython/profile_default/文件夹中，这里default是profile 的用户名，也是我们打开notebook的默认用户。但是有时候我们会发现这里面的东西并不像配置文件。其实我们只需要输入如下命令即可把默认的配置文件给显示出来：
$ipython profile create default
很好理解，就是显式的创建一个default用户（当然我们还可以创建其他名字的用户）。

现在，打开~/.ipython/profile_default/ipython_notebook_config.py文件，这个就是我们的配置文件了。

设置密码

由于ipython notebook 远程登陆的话权限很大，因此还是设置一个密码来保证安全。

1、首先我们要得到密码的哈希值，在ipython 中进行如下操作:

1
2
3
4
5
In [1]: from IPython.lib import passwd
In [2]: passwd()
Enter password:
Verify password:
Out[2]: 'sha1:7467b7351f12:79fc6599f863a8daf048c655f60adf03003b87f3'

在这里输入密码，就可以得到他的hash值了。记住Out[2]中的字符串。

2、然后在配置文件ipython_notebook_config.py中写如入下面的信息：

1
2
c = get_config()
c.NotebookApp.password =u'sha1:7467b7351f12:79fc6599f863a8daf048c655f60adf03003b87f3'

这样就可以了。

配置接入点

这是稍微麻烦的地方，因为ipython的端口默认是8888而且不能跟apache2的冲突，所以在配置的时候需要想办法用一个优雅的姿势来用apache2访问到8888端口。

1、首先对配置文件进行一下设置，打开ipython_notebook_config.py，加入下面的信息：

1
2
3
4
5
6
7
8
9
10
11
c.IPKernelApp.pylab = 'inline'  # if you want plotting support always
c.NotebookApp.ip = '*'
c.NotebookApp.open_browser = False
# It is a good idea to put it on a known, fixed port
c.NotebookApp.port = 8888
````
很明显，我们可以自由的配置端口。由于ipython notebook会默认打开一个浏览器窗口，而服务器是打不开浏览器的，所以我么会把c.NotebookApp.open_browser设置为false防止warning。

c.NotebookApp.ip经检测还是设置为*好。。。否则会出现很多莫名的错误。

2、为了让apache2能自然的访问到8888端口，我们配置下.htaccess文件。打开网站根目录下的.htaccess文件，添加下面的语句：

RewriteCond %{HTTP_HOST} ^ipython.mythsman.com$ [NC]
RewriteCond %{REQUEST_URI} ^/(.)$ [NC]
RewriteRule ^(.)$ http://ipython.mythsman.com:8888/%1 [R=301,L]

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
意思很简单，就是把ipython.mythsman.com的80端口重定向到8888端口。而这个站点的根目录其实就是我们打开ipython notebook 服务的那个目录。

事实上，我们可以在当前站点的任何子域名下输入8888端口进入这里。

## Jupyter

不知道为什么ipython notebook对于文件的同步性并不好。相比之下Jupyter做的就好多了，而且界面也更加漂亮。其实我个人觉得Jupyter就是ipython notebook的一种更新。

详细信息参见[官方文档](http://jupyter.readthedocs.org/en/latest/)。

### 安装

用`$pip install jupyter`等一系列类似手段，不解释了。



### 配置

#### 配置文件

其实jupyter本身的配置还是有点麻烦的，不过我们可以把ipython notebook的配置文件移植过来～～参考下面的映射表：

|ipython notebook | jupyter|
|-|-|
|~/.ipython/profile_default/ipython_notebook_config.py|~/.jupyter/jupyter_notebook_config.py|
|~/.ipython/profile_default/static/custom|~/.jupyter/custom|
|~/.ipython/profile_default/ipython_nbconvert_config.py|~/.jupyter/jupyter_nbconvert_config.py|
|~/.ipython/profile_default/ipython_qtconsole_config.py|~/.jupyter/jupyter_qtconsole_config.py|
|~/.ipython/profile_default/ipython_console_config.py|~/.jupyter/jupyter_console_config.py|

主要是第一个配置文件了，参照ipython的配置文件的写法照着搬过来就可以了，懒得去看他的配置命令。。


#### 配置服务器

与ipython一样配置下.htaccess即可。

RewriteCond %{HTTP_HOST} ^jupyter.mythsman.com$ [NC]
RewriteCond %{REQUEST_URI} ^/(.)$ [NC]
RewriteRule ^(.)$ http://jupyter.mythsman.com:8888/%1 [R=301,L]

在想要显示的文件夹处输入jupyter notebook，即可在jupyter.mythsman.com处访问到jupyter。

Posted on | In

由于有关caffe的开发资料实在太少，单单是这个问题就困扰了我半天。最后终于找到了一个大腿—-beenfrog，也是一个正在学习caffe框架的研究人员。博客mark下，以后有问题可以去这里找。

本文代码部分主要参考于此文。

简述

lmdb是一种轻量级的数据库，caffe中主要就是使用lmdb模块来进行图像数据集的保存。据说是因为lmdb有读取速度快，支持多线程、多进程并发，等这样那样的优点（具体见官网，虽然我暂时没有看出来，据我所知网上查找lmdb文档的人大都仅仅是为了使用caffe的），注意到这个数据库其实并没有任何压缩处理的作用，他的目的只是为了快速的索引和存取。他的数据都会带着一定的数据结构从而使的体积略微增大。

事实上如果仅仅看lmdb的用法是无法直接应用于图像文件的处理的。由于caffe是将图像以他自带的数据类型的形式传入lmdb中的，因此我们必须结合caffe的数据类型才能完成读取和使用。

lmdb的安装

参考官方文档，安装如下依赖即可：

1
2
3
4
udo apt-get install cython
sudo apt-get install libffi-dev python-dev build-essential
sudo apt-get install python-cffi
sudo easy_install lmdb

实践再次映证了用easy_install 安装要比pip好的多(在版本跟的上的情况下)。

生成数据文件

这里以zhengfang/文件夹下的1.png至1000.png为例

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
#coding:utf-8
import lmdb
import numpy as np
import cv2
import caffe
from caffe.proto import caffe_pb2

#basic setting
lmdb_file = 'lmdb_data'#期望生成的数据文件
batch_size = 200       #lmdb对于数据进行的是先缓存后一次性写入从而提高效率，因此定义一个batch_size控制每次写入的量。

# create the leveldb file
lmdb_env = lmdb.open(lmdb_file, map_size=int(1e12))#生成一个数据文件，定义最大空间
lmdb_txn = lmdb_env.begin(write=True)              #打开数据库的句柄
datum = caffe_pb2.Datum()                          #这是caffe中定义数据的重要类型

for x in range(1000):
    x+=1
    img=cv2.imread('zhengfang/'+str(x)+'.png').convert('RGB')     #从zhengfang/文件夹中依次读取图像

    # save in datum
    data = img.astype('int').transpose(2,0,1)      #图像矩阵，注意需要调节维度
    #data = np.array([img.convert('L').astype('int')]) #或者这样增加维度
    label = x                                      #图像的标签，为了方便存储，这个必须是整数。
    datum = caffe.io.array_to_datum(data, label)   #将数据以及标签整合为一个数据项

    keystr = '{:0>8d}'.format(x-1)                 #lmdb的每一个数据都是由键值对构成的，因此生成一个用递增顺序排列的定长唯一的key
    lmdb_txn.put( keystr, datum.SerializeToString())#调用句柄，写入内存。

    # write batch
    if x % batch_size == 0:                        #每当累计到一定的数据量，便用commit方法写入硬盘。
        lmdb_txn.commit()
        lmdb_txn = lmdb_env.begin(write=True)      #commit之后，之前的txn就不能用了，必须重新开一个。
        print 'batch {} writen'.format(x)

lmdb_env.close()                                   #结束后记住释放资源，否则下次用的时候打不开。。。

输出：

1
2
3
4
5
batch 200 writen
batch 400 writen
batch 600 writen
batch 800 writen
batch 1000 writen

照着注释，很好理解了。值得一提的是，生成的文件是以文件夹的形式保存的，内部由data.mdb lock.mdb组成，不用管他。

这里需要强烈注意的一点是，datum里的图像是需要有长宽高三个维度的，而且第一个维度默认是通道数，而我们平常的图像格式是第三位是通道。因此我们需要在读入普通图像的时候将他转置一下存入datum中，或者在外面增加一层。

读取数据文件

读取上面生成的数据文件。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
#coding:utf-8

import caffe
import lmdb
import numpy as np
import cv2
from caffe.proto import caffe_pb2

lmdb_env = lmdb.open('lmdb_data')#打开数据文件
lmdb_txn = lmdb_env.begin()      #生成句柄
lmdb_cursor = lmdb_txn.cursor()  #生成迭代器指针
datum = caffe_pb2.Datum()        #caffe定义的数据类型

for key, value in lmdb_cursor:   #循环获取数据
    datum.ParseFromString(value) #从value中读取datum数据

    label = datum.label          #获取标签以及图像数据
    data = caffe.io.datum_to_array(datum)
    print data.shape
    print datum.channels
    image =data.transpose(1,2,0)
    cv2.imshow('cv2.png', image) #显示
    cv2.waitKey(0)

cv2.destroyAllWindows()
lmdb_env.close()

输出：（图像略）

1
2
(3, 27, 72)
3

基本是和写入程序一一相对，很好理解。同样需要注意的是如需要显示图像，则需要将数据转置回来。

Posted on | In

从python2.6之后，python中的字符串就有了str.format()函数这一格式控制的强大工具。相比于之前使用%的格式控制手段，str.format()函数显然更加符合我们的思维习惯，而且更加简洁。

语法

作为字符串的一个方法，它以{}和:来代替%，进行格式控制。

定位

通过位置

1
2
3
4
5
6
In [1]: '{0},{1}'.format('kzc',18)
Out[1]: 'kzc,18'
In [2]: '{},{}'.format('kzc',18)
Out[2]: 'kzc,18'
In [3]: '{1},{0},{1}'.format('kzc',18)
Out[3]: '18,kzc,18'

十分好理解，就是在字符串中用{}中的值来指定format中用以代替他的值。可以接受不限个参数，位置可以不按顺序，可以不用或者用多次。

通过关键字参数

1
2
In [5]: '{name},{age}'.format(age=18,name='kzc')
Out[5]: 'kzc,18'

就是可以用键值对的形式给参数列表中的元素赋值。

通过对象属性

1
2
3
4
5
6
7
8
class Person:
  def __init__(self,name,age):
    self.name,self.age = name,age
    def __str__(self):
      return 'This guy is {self.name},is {self.age} old'.format(self=self)

In [2]: str(Person('kzc',18))
Out[2]: 'This guy is kzc,is 18 old'

在参数列表中指定要显示的对象的属性，并在format中传入这个对象。

通过下标

1
2
3
In [7]: p=['kzc',18]
In [8]: '{0[0]},{0[1]}'.format(p)
Out[8]: 'kzc,18'

结合位置与下表进行传参。

格式限定符

填充与对齐

填充常跟对齐一起使用
^、<、>分别是居中、左对齐、右对齐，后面带宽度
:号后面带填充的字符，只能是一个字符，不指定的话默认是用空格填充
比如

1
2
3
4
5
6
7
In [15]: '{:>8}'.format('189')
Out[15]: '   189'
In [16]: '{:0>8}'.format('189')
Out[16]: '00000189'
In [17]: '{:a>8}'.format('189')
Out[17]: 'aaaaa189'

精度与类型

1
2
In [44]: '{:.2f}'.format(321.33345)
Out[44]: '321.33'

精度常跟类型f一起使用，本例中.2表示长度为2的精度，f表示float类型。

进制转换

1
2
3
4
5
6
7
8
In [54]: '{:b}'.format(17)
Out[54]: '10001'
In [55]: '{:d}'.format(17)
Out[55]: '17'
In [56]: '{:o}'.format(17)
Out[56]: '21'
In [57]: '{:x}'.format(17)
Out[57]: '11'

b、d、o、x分别是二进制、十进制、八进制、十六进制，跟在冒号后面

分割字符

1
2
In [47]: '{:,}'.format(1234567890)
Out[47]: '1,234,567,890'

用逗号还能用来做金额的千位分隔符。

Posted on | In

在整理最近学习的知识时，突然意识到机器学习与深度学习的差别。之前学的gd算法，logistic回归算法，svm算法等都属于机器学习的范畴，而深度学习与他们其实是并列的关系，同属机器学习这个大的范畴。个人认为深度学习其实是曾经没落的BP神经网络的发展，可以说神经网络就是趁着深度学习的浪潮借尸还魂的。不过，深度学习想要理解起来也不是那么容易的。既然如此，我们何不搭建一个黑盒环境，先直观感受一下深度学习的魅力再慢慢研究呢？这里我们用到的工具就是Caffe深度学习框架。

简介

说道Caffe，我们不得不提他的作者贾扬清（点击进入个人主页）。没错是一个中国人，而且是一个非常厉害的中国人（我词穷）。话不多说，我们来看看他的CV：

Yangqingcv

只是很少的一部分，不说话，我们用心来感受就可以了。

Caffe其实是他与他的朋友利用课余时间写的一个框架，而且写的时候正逢他写博士毕业论文的时候。这么推来大概是2014年左右完成的。

因此说白了Caffe只是一个普通的小团体的开源项目而已，就连他的官网也是挂靠在berkeleyvision.org上。官网：http://caffe.berkeleyvision.org/

环境搭建

根据惯例，本人还是一切以官方文档做参考，搭建在ubuntu14.04系统上。

下载

项目托管在github上，见他的github即可。

依赖

依赖有很多，除了常见的包之外，还要下这些：

1
2
3
4
sudo apt-get install libprotobuf-dev libleveldb-dev libsnappy-dev libopencv-dev libhdf5-serial-dev protobuf-compiler
sudo apt-get install --no-install-recommends libboost-all-dev
sudo apt-get install libgflags-dev libgoogle-glog-dev liblmdb-dev
sudo apt-get install libatlas-base-dev

当然，原则上还要有GPU的依赖，但是不用GPU照样可以工作，只是要注意之后的参数设定。

编译

和OpenCV一样，这里用cmake编译。

打开项目根目录，依次输入以下指令即可：

1
2
3
4
5
6
mkdir build
cd build
cmake ..
make all
make install
make runtest

这样就会在build文件加里面生成很多的文件。注意，这些文件不能删的，因为事实上他编译得到的头文件啊、二进制文件啊都在这里，将来使用caffe的时候都要用到。

如果需要安装python接口的话还要输入make pycaffe，而且还需要下载python的依赖。在项目根目录下，输入：

1
2
cd ./build/install/python/
for req in $(cat requirements.txt); do sudo pip install $req; done

就是要满足这个文件里的所有依赖(除此之外还有一个skimage.so需要下载i)。但是，由于国内的pip镜像普遍不全，如果不翻墙，总有一些包下不下来。因此，我手动打开requirements.txt，查看具体需要的包，然后通过$sudo apt-get install python-*****的办法下载某些无法下载的包。

如果完全下载好了，那么在项目中进入./python目录中打开ipython，就能import caffe 了。当然，为了方便，我们可以把caffe文件夹（caffe模块的位置）复制到/usr/lib/python2.7/dist-packages 目录下，这样就可以在任意的地方调用了。

运行测试

事实上在编译阶段就进行了测试，只是这样并没有具体的感悟，下面就结合mnist来简单测试下。

在项目根目录下输入：

1
2
./data/mnist/get_mnist.sh
./examples/mnist/create_mnist.sh

表示下载mnist训练集，并转换为指定格式。

然后需要打开./examples/mnist/lenet_solver.prototxt 文件，并将最后一行的GPU改成CPU。(毕竟之前没有装GPU，果断是不支持的)。

最后就可以跑训练程序了。

1
./examples/mnist/train_lenet.sh

得到结果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
......
I0227 16:58:31.228276  5528 solver.cpp:228] Iteration 9800, loss = 0.0148246
I0227 16:58:31.228313  5528 solver.cpp:244]     Train net output #0: loss = 0.0148246 (* 1 = 0.0148246 loss)
I0227 16:58:31.228335  5528 sgd_solver.cpp:106] Iteration 9800, lr = 0.00599102
I0227 16:58:37.839028  5528 solver.cpp:228] Iteration 9900, loss = 0.00598603
I0227 16:58:37.839066  5528 solver.cpp:244]     Train net output #0: loss = 0.00598603 (* 1 = 0.00598603 loss)
I0227 16:58:37.839088  5528 sgd_solver.cpp:106] Iteration 9900, lr = 0.00596843
I0227 16:58:44.390777  5528 solver.cpp:454] Snapshotting to binary proto file examples/mnist/lenet_iter_10000.caffemodel
I0227 16:58:44.395431  5528 sgd_solver.cpp:273] Snapshotting solver state to binary proto file examples/mnist/lenet_iter_10000.solverstate
I0227 16:58:44.425786  5528 solver.cpp:317] Iteration 10000, loss = 0.00296525
I0227 16:58:44.425813  5528 solver.cpp:337] Iteration 10000, Testing net (#0)
I0227 16:58:48.531802  5528 solver.cpp:404]     Test net output #0: accuracy = 0.9909
I0227 16:58:48.531833  5528 solver.cpp:404]     Test net output #1: loss = 0.0281938 (* 1 = 0.0281938 loss)
I0227 16:58:48.531841  5528 solver.cpp:322] Optimization Done.
I0227 16:58:48.531859  5528 caffe.cpp:222] Optimization Done.

最后可以发现拟合的准确率为99.09%。

至于每一个命令文件到底是做了什么可以直接打开看看。

Posted on | In

在进行图像检测或者是识别的时候，我们需要提取出一些有特征的点加以识别，最常用的就是基于点的识别。这里所谓的点，其实就是一些重要的点，比如轮廓的拐角，线段的末端等。这些特征比较容易识别，而且不容易受到光照等环境的影响，因此在许多的特征匹配算法中十分常见。

常见的特征点提取算法有Harris算 子（改进后的Shi-Tomasi算法）、Moravec算子、Forstner算子、小波变换算子等。现在就先介绍一下最常用的Harris角点检测算法。

简介

Harris算法的思想很简单，也很容易理解。

我们知道角点附近的区域相比于其他地方有这样一个显著的特点，就是无论沿着哪一个方向看，他灰度的变化率始终是很大的。也就是说假设我们有一个矩形窗口罩在角点附近，将这个窗口顺着任意方向移动一小段距离得到一个新的区域，将这个新的区域与旧的区域对应点的灰度做差得到的值始终很大。相比之下，平滑区域的变化就很小，而边缘区域沿着某些方向变化率大、某些方向变化率小。

Harris算法利用的就是这个特点，他首先定义了一个窗口函数$w(x,y)$来表示他选择的窗口区域，$(x,y)$表示点的坐标，$w(x,y)$表示这个坐标所占的权值。有时候我们用0-1赋值，表示选定一块区域，也有时候我们用高斯滤波减少噪点的影响；然后定义了一个方向向量$(u,v)$，以及$E(u,v)$表示窗口沿着$(u,v)$方向移动后的梯度变化情况。并将$E(u,v)$做如下定义：

$E(u,v)=\underset{x,y}{\Sigma} w(x,y)[I(x+u,y+v)-I(x,y)]^2$

其中$I(x,y)$表示点$(x,y)$的灰度值。

根据上面的介绍我们知道角点的特征就是$E(u,v)$的值取较大值。那么为了更方便的计算，我们对他需要进行一下化简：

泰勒展开：$I(x+u,y+v)-I(x,y)=I(x,y)+uI_x+vI_y$

其中$I_x,I_y$分别为灰度沿x,y方向的导数。

转化为矩阵形式：$E(u,v)=[u\ v]M\begin{bmatrix}u\\v\end{bmatrix},M=\underset{x,y}\Sigma w(x,y)\begin{bmatrix}I_x^2&I_xI_y\\I_xI_y&I_y^2\end{bmatrix}$

最后定义一个估价函数R：

$R=det(M)=k(trace(M))^2$

其中$det(M)=\lambda_1\lambda_2,trace(M)=\lambda_1+\lambda_2$,k是一个控制参数。

$\lambda_1,\lambda_2$为M的特征值。

这个估价函数个特性，就是当R较小时，图像是平坦的；当R小于０时，图像是一个边缘；当R很大时，这个图像是一个角点。因此通常我们会对R设置一个阈值，大于这个阈值的点我们可以看做是角点。

OpenCV调用

OpenCV里封装了Harris算法，调用下看看效果，就不自己实现了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#coding:utf-8
import cv2
import numpy as np

img=cv2.imread('test.png',cv2.IMREAD_GRAYSCALE)
cv2.imshow('gray.png',img)#控制背景为黑色
cv2.imwrite('gray.png',img)
im=np.float32(img)
dst=cv2.cornerHarris(im,3,3,0.04)#生成估价矩阵

height,width=dst.shape
mark=dst>0.01*dst.max()#寻找具有较大权值的像素点
height,width=mark.shape
for i in xrange(height):
    for j in xrange(width):
        if mark[i][j]:
            cv2.circle(img,(j,i),5,255)#标记得到的点
cv2.imshow('Harris.png',img)
cv2.imwrite('Harris.png',img)
cv2.waitKey(0)
cv2.destroyAllWindows()

重要的步骤就是在估价矩阵里找到较大的那些点。

效果图

gray
Harris

Posted on | In

击中击不中变换（Hit Miss Transform ，HMT），是通过同时探测图像的内部和外部，进而获取更多的内外标记，体现更多信息的一个方法。他的应用有很多，特别是在图像识别以及图像细化方面。

定义

既然既要有击中也要有击不中，那么显然我们需要两个结构基元E和F，我们把这两个结构基元记为一个结构元素对$B=(E,F$。其中一个用来探测图像内部，一个用来探测图像外部。

对于给定的图像A，我们定义用B对他进行的击中击不中变换为

$A*B=(A\Theta E)\bigcap (A^c\Theta F)$

即用E对A进行腐蚀，用F对A的补集进行腐蚀，并将得到的结果求交集。

这样说可能冠冕堂皇了一点，简单的讲其实就在判断图像中的某一个像素，将这个像素与结构基元的‘原点’对应，看看结构基元中的‘击中’基元是否能完全被图像覆盖、结构基元中的‘击不中部分’是否能与图像没有任何交集。如果答案均为是，那么这个点就可以得到保留，否则就舍弃。当然，这里用了腐蚀的方法使得图像更加容易被识别进去。

说白了一点其实就是相当于在用一个模子来测试这个图片，并将所有符合这个模子的中心点标记下而已。

应用

物体识别

这个很显而易见了，只要将击中基元设计成为需要识别的图像的样子，击不中基元设计为击中基元取反后的样子即可。然后用这两个基元对图像进行检测即可。

不过有时候图像会发生略微变形，为了提高识别率，我们通常会将两个结构基元都稍微腐蚀一下，给图像变形留下些余地。

思路比较清楚，具体代码依据实际情况编写即可，这里不做举例。

细化

定义

图像细化是一种最常见的使用击中击不中变换的形态学算法。他的目的就如他的字面意思，将图像变细，就像取到了图像的骨架。

其定义是，对于图像S，利用结构对$B=(E,f)$进行细化记为：

$S\bigotimes B=S\setminus (S*B)$

其中$S*B$表示用B对S进行细化得到的标记点，$\setminus$表示求差集。

下面介绍一种细化的方法。

方法

设置八个分别代表不同边界的击中击不中模板：

$\begin{bmatrix}\bullet&\bullet&\bullet\\&\circ&\\\circ&\circ&\circ\end{bmatrix}$$\begin{bmatrix}&\bullet&\bullet\\\circ&\circ&\bullet\\\circ&\circ&\end{bmatrix}$$\begin{bmatrix}\circ&&\bullet\\\circ&\circ&\bullet\\\circ&&\bullet\end{bmatrix}$$\begin{bmatrix}\circ&\circ&\\\circ&\circ&\bullet\\&\bullet&\bullet\end{bmatrix}$

$\begin{bmatrix}\circ&\circ&\circ\\&\circ&\\\bullet&\bullet&\bullet\end{bmatrix}$$\begin{bmatrix}&\circ&\circ\\\bullet&\circ&\circ\\\bullet&\bullet&\end{bmatrix}$$\begin{bmatrix}\bullet&&\circ\\\bullet&\circ&\circ\\\bullet&&\circ\end{bmatrix}$$\begin{bmatrix}\bullet&\bullet&\\\bullet&\circ&\circ\\&\circ&\circ\end{bmatrix}$

其中$\circ$表示击中模板，$\bullet$表示击不中模板，０表示不设置限制。每个矩阵的中心位置设为原点。

接下来我们要做的就是用着八个模板依次对待处理的图像做击中击不中变换，每轮操作都会使图像变细。经过多次的迭代之后图像保持不变，这个结果就是经过细化后的图像了。

由于OpenCV中没有直接进行HMT的模板，因此手敲了一个（大概意思差不多是这样，就是效率不敢恭维。。）

OpenCV实现

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
#coding:utf-8
import cv2
import numpy as np

im=cv2.imread('test.png',cv2.IMREAD_GRAYSCALE)
thresh,im=cv2.threshold(im,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
cv2.imshow('binary.png',im)#控制背景为黑色
cv2.imwrite('binary.png',im)

arr=np.array(im)
mat=[]#HMT模板矩阵
mat.append(np.array([[-1,-1,-1],[0,1,0],[1,1,1]]))
mat.append(np.array([[0,-1,-1],[1,1,-1],[1,1,0]]))
mat.append(np.array([[1,0,-1],[1,1,-1],[1,0,-1]]))
mat.append(np.array([[-1,-1,-1],[0,1,0],[1,1,1]]))
mat.append(np.array([[1,1,0],[1,1,-1],[0,-1,-1]]))
mat.append(np.array([[1,1,1],[0,1,0],[-1,-1,-1]]))
mat.append(np.array([[-1,0,-1],[-1,1,1],[-1,0,1]]))
mat.append(np.array([[-1,-1,0],[-1,1,1],[0,1,1]]))

height,width=arr.shape

while True:#迭代直至无变化
    before=arr.copy()
    for m in mat:#使用八个模板进行变换
        mark=[]
        for i in xrange(height-2):#对每个非边界点进行测试
            for j in xrange(width-2):
                reg=True
                for im in xrange(3):
                    for jm in xrange(3):
                        if not arr[i+1][j+1]==255:
                            continue
                        if m[im][jm]==1 and arr[i+im][j+jm]==0:
                            reg=False
                        if m[im][jm]==-1 and arr[i+im][j+jm]==255:
                            reg=False
                if reg:#找到标记，删除
                    mark.append((i+1,j+1))
        for it in mark:
            x,y=it
            arr[x][y]=0
    if (before==arr).all():
        break

cv2.imshow('thin.png',arr)
cv2.imwrite('thin.png',arr)
cv2.waitKey(0)
cv2.destroyAllWindows()

就是用八个模板不停的求卷积，并将标记到的点删除，直至无法删除为止。

运行结果

binary

thin

Posted on | In

开闭运算是二值形态学的重要部分，是对腐蚀和膨胀算法的扩展应用，在图像的去噪方面也是十分的常用。

二值开运算

定义

用结构元素B（即模板核）对图像A进行开运算，可以用符号$A\circ B$来表示，其定义式为：

$A\circ B=(A\Theta B)\bigoplus B$

所以开运算实际上就是A先被B腐蚀，再被B膨胀。

作用

很明显，开运算可以消除一些很小的背景噪点（椒盐噪声），平滑较大物体的边界而不明显改变其体积。也会磨光矩形的内边缘，以及分离一些粘连目标。

（这里所说的背景噪点是指以黑色为背景色，白色为前景色）

OpenCV调用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#coding:utf-8
import cv2
import numpy as np

im=cv2.imread("test.png",cv2.IMREAD_GRAYSCALE)#以二值图读入

im=cv2.bitwise_not(im)#由于背景为白色，将其置反从而统一标准

thresh,im=cv2.threshold(im,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)#OTSU转化为二值图

cv2.imshow('binary',im)
cv2.imwrite('binary.png',im)

kernel=np.ones((2,2),np.uint8)#模板核

opening=cv2.morphologyEx(im,cv2.MORPH_OPEN,kernel)#进行开运算消除背景噪声

cv2.imshow('opening',opening)
cv2.imwrite('opening.png',opening)

cv2.waitKey(0)
cv2.destroyAllWindows()

这里调用的是cv2.morphologyEx函数，这是专门用于处理形态学问题的，传入适当参数即可。

由于图像本身是以白色为背景，因此将他的像素取个反，使他符合约定俗称的黑色背景。

（写的时候发现一个问题，就是imread好像并不直接支持gif图像格式的读取，因此本来是test.gif的图片还得事先进行convert，转化为test.png然后才好使用）

效果图

binary

opening

二值闭运算

定义

用结构元素B（即模板核）对图像A进行闭运算，可以用符号$A\cdot B$来表示，其定义式为：

$A\cdot B=[A\bigoplus (-B)]\Theta (-B)$

所以闭运算实际上就是A先被B膨胀，再被B腐蚀。

作用

闭运算在去除图像的前景噪声方面具有很好的应用，通过闭运算之后，图像原有目标的间断可以得到连接、目标内部的孔洞可以得到填充而基本不改变原图的大小以及形态。同时也会磨光凸向图像内部的边角。

OpenCV调用

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#coding:utf-8
import cv2
import numpy as np

im=cv2.imread("test.png",cv2.IMREAD_GRAYSCALE)#以二值图读入

thresh,im=cv2.threshold(im,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)#OTSU转化为二值图

cv2.imshow('binary',im)
cv2.imwrite('binary.png',im)

kernel=np.ones((2,2),np.uint8)#模板核

closing=cv2.morphologyEx(im,cv2.MORPH_CLOSE,kernel)#进行闭运算消除前景噪声

cv2.imshow('closing',closing)
cv2.imwrite('closing.png',closing)

cv2.waitKey(0)
cv2.destroyAllWindows()

没啥好说的了，对着开操作改下参数就好了。

效果图

binary
closing

Posted on | In

数学形态学（Mathematical morphology） 是一门建立在格论和拓扑学基础之上的图像分析学科，是数学形态学图像处理的基本理论。其基本的运算包括：二值腐蚀和膨胀、二值开闭运算、骨架抽取、极限腐蚀、击中击不中变换、形态学梯度、Top-hat变换、颗粒分析、流域变换、灰值腐蚀和膨胀、灰值开闭运算、灰值形态学梯度等。

这是我们在进行图像处理时十分常见的算法工具，在实际中有着非常广泛的应用。

下面就先介绍下二值形态学中的两个重要操作－－腐蚀与膨胀。

二值腐蚀

简述

所谓二值腐蚀，就是对一个二值图进行腐蚀操作。我们对一个二值图像进行腐蚀操作首先需要一个模板核。模板核由一个二值矩阵组成，当中还定义了一个原点，表示核的核心。

接下来，对于二值图像的每一个像素点（无论０还是１），我们将他与模板核的核心对齐，查看模板核中所有值为１的格点组成的图像是否被完全包含在那个二值图像中。如果是，则将该像素点更新为１，否则，更新为０。

通俗的讲，其实就是用我们的模板来对图像进行检测，以便找出图像内部可以放得下该基元的区域。

这个过程的效果是消除边界点，使得边界向内部收缩。

定义

集合$A$被$B$腐蚀，表示为$A\Theta B={c|B+c\subset A }$

其中$A$为待腐蚀的图像，$B$为模板向量。

比如二值图像：

$S=\begin{bmatrix}0&1&0&1&0\\0&1&1&0&1\\0&1&1&1&0\end{bmatrix}$

以及模板：

$E=\begin{bmatrix}1&0\\1&1_\triangle\end{bmatrix}$

其中$\ _\triangle$表示模板的原点。

那么$E$对$S$腐蚀的结果就为：

$S\Theta E=\begin{bmatrix}0&0&0&0&0\\0&0&1&0&0\\0&0&1&1&0\end{bmatrix}$

OpenCV调用

本来想用pil实现下的，然而突然想通了，既然OpenCV有实现好的函数，为何还要用pil折腾呢？于是就用OpenCV的python接口体验下效果：(python接口的安装只需直接apt-get install python-opencv即可)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#coding:utf-8
import numpy as np
import cv2

im=cv2.imread('test.jpg',cv2.IMREAD_GRAYSCALE)#读取图像

thresh,im=cv2.threshold(im,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)#用OTSU自动获取阈值并进行二值化，第一个参数表示图像，第二个表示设置阈值（由于我们用OTSU自动设置，所以这里必须填０），第三个参数表示将超过正常范围的像素设置的值，最后一个传入控制参数。

cv2.imshow('binary',im)
cv2.imwrite('binary.png',im)

kernel=np.ones((6,6),np.uint8)#生成一个6x6的核

erosion=cv2.erode(im,kernel,iterations=1)#调用腐蚀算法

cv2.imshow('erosion',erosion)
cv2.imwrite('erosion.png',erosion)
cv2.waitKey(0)
cv2.destroyAllWindows()

图样

binary

erosion

二值膨胀

简述

晓得了二值腐蚀，那么二值膨胀也就很好理解了。他的作用和腐蚀恰好相反，他能够用来填补一些小洞洞，将两个物体连接起来。

最终效果就是扩展边界，使得图像变得膨胀。对于膨胀，其实可以用腐蚀来理解。对图像的膨胀其实也就是对背景的腐蚀。

当然也可以用类似腐蚀的方法来描述，我们需要用一个定义了原点的模板核，这个核是一个二值矩阵。然后对于待处理图像的每一个值为１的像素点，将他与模板核的原点对其，然后对于模板核中值为１的点在图像中对应的像素更新为１。

通俗的讲，就是将每一个像素点用模板核来代替，并将所得到的所有图像求一个并操作。

定义

集合$A$被模板$B$膨胀，表示为$A\bigoplus B=[A^c\Theta (-B)]^c$

其中$A^c$表示A的补集，$-B$表示将$B$旋转180°。

以上是用腐蚀进行的定义，还可以直接定义：

$A\bigoplus B=\cup{A+b|b\in B}$

这个式子也被称作明夫斯基和形式。

OpenCV调用

同之前一样，调用下OpenCV的接口体验下效果：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#coding:utf-8
import numpy as np
import cv2
im=cv2.imread('test.jpg',cv2.IMREAD_GRAYSCALE)#读取图像

thresh,im=cv2.threshold(im,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)#用OTSU自动获取阈值并进行二值化，第一个参数表示图像，第二个表示设置阈值（由于我们用OTSU自动设置，所以这里必须填０），第三个参数表示将超过正常范围的像素设置的值，最后一个传入控制参数。

cv2.imshow('binary',im)
cv2.imwrite('binary.png',im)

kernel=np.ones((6,6),np.uint8)#生成一个6x6的核

dilation=cv2.dilate(im,kernel,iterations=1)#调用膨胀算法

cv2.imshow('dilation',dilation)
cv2.imwrite('dilation.png',dilation)
cv2.waitKey(0)
cv2.destroyAllWindows()

图样

binary

dilation

Posted on | In

OTSU算法又叫最大类间方差阈值分割算法，也叫大津算法，是在1980年由日本的大津展之提出，是由最小二乘法推导而来，用于一些简单的阈值确定。

对于一个灰度图，我们有时候非常想把他用一个阈值将他的前景和背景区分开来。我们可以合理的假设为如果将图像的像素分布图画出来，那么图像上应该有两个峰，即前景色和背景色。在这两个峰之间肯定有一个谷，那么我们就可以将阈值设在这里，从而对图像达到一个良好的分割效果。

怎样确定这个阈值呢？OTSU算法说，我们可以求出用这个阈值分割后的两个图像的类间方差。对于每一个可能的阈值，我们计算并取出类间方差最大的那个像素值，此时这个值就可以较好的对图像进行分割。

算法

１、将灰度值分为$0-m$，对于$0-m$的每一个灰度$t$，将他作为阈值将图像分割为灰度为$0-t$以及$t+1-m$这两部分。

２、计算每一部分的所占比例$w_1,w_2$，每一部分的平均灰度值$u_1,u_2$，以及总的平均灰度值$u$。

３、计算他们的类间方差$\delta^2=w_1(u_1-u)^2+w_2(u_2-u)^2=w_1w_2(u_1-u_2)^2$

４、取出类间方差最大时对应的阈值ｔ，这就可以作为我们最终所取的阈值。

测试

用python写的简单实现：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
import Image
import numpy as np

im=Image.open('test.png')
im.show()
im=im.convert('L')
im.show()
arr=np.array(im)
width,height=arr.shape
arr=arr.reshape(width*height)

freq=np.zeros(256)*1.0

total=0.
point=0.
for i in arr:
    freq[i]+=1
    total+=i
    point+=1

u=total/point

w1=0.
w2=1.
u1=0.
u2=u
eps=0.
threshold=0

for i in range(255):
    if freq[i]==0 or w2*point-freq[i]==0:
        continue
    u1=(u1*w1*point+i*freq[i])/(w1*point+freq[i])
    u2=(u2*w2*point-i*freq[i])/(w2*point-freq[i])
    w1=w1+freq[i]/point
    w2=w2-freq[i]/point
    eps_now=w1*(u1-u)*(u1-u)+w2*(u-u2)*(u-u2)
    if(eps_now>eps):
        eps=eps_now
        threshold=i

table=[]
for i in range(256):
    if i>threshold :
        table.append(255)
    else:
        table.append(0)
im=im.point(table,'L')
im.show()

效果图：