Resive world

Come world to Record life


  • Home

  • Tags

  • Categories

  • Archives

  • Sitemap

  • Search

最小二乘法

Posted on 2016-02-22 | In Machine Learning

最小二乘法,说白了其实就是解决线性回归问题的一个算法。这个算法最早是由高斯和勒让德分别独立发现的,也是当今十分常见的线性拟合算法,并不复杂。

我们常用的最小二乘法有两种,一种是普通方程表示的简单线性拟合问题,另一种是矩阵表示的高维度的线性拟合问题。

普通最小二乘法

他解决的基本问题其实就是给定一些数对$(x_1,y_1),(x_2,y_2),(x_3,y_3),….$,让你求出参数$\beta_0,\beta_1$,使得直线$y=\beta_0+\beta_1 x$能够最好的拟合这个数据集,也就是使得他的平方损失函数取到最小值,即
$$Q=\underset{i=1}{\overset{n}{\sum}} (y_i-\beta_0-\beta_1x_i)^2$$
最小。

为了取得最小值,我们分别对$\beta_0,\beta_1$求偏导,并使之为0。

$$
\left\{\begin{aligned}&\frac{\partial Q}{\partial \beta_0}=-2\underset{i=1}{\overset{n}{\sum}}(y_i-\beta_0-\beta_1 x_i)=0\\& \frac{\partial Q}{\partial \beta_1}=-2\underset{i=1}{\overset{n}{\sum}}x_i(y_i-\beta_0-\beta_1 x_i)=0\end{aligned}\right.
$$

解得:

$$\left\{\begin{aligned}&\beta_0=\frac{n\sum x_iy_i-\sum x_i\sum y_i}{n\sum x_i^2-(\sum x_i)^2}\\& \beta_1=\frac{\sum x_i^2\sum y_i-\sum x_i\sum x_iy_i}{n\sum x_i^2-(\sum x_i)^2}\end{aligned}\right.$$

套用这个公式得到的参数$\beta_0,\beta_1$就是最好的拟合参数了。

矩阵最小二乘法

用矩阵表示的最小二乘法则更加方便,能够用非常简单的矩阵形式进行计算,而且能拟合多维度的线性方程。

对于线性回归,我们要做的事情其实可以近似等同于解线性方程AX=Y,其中A是mn的矩阵,X,Y是1m的矩阵。m是数据的对数,n是数据的维数加1(因为还有常数),而且n应该小于m。

当然,这个方程是没有解得,不过我们可以通过数值运算计算出一个解,可以证明这个解就是平方损失函数最小的解。

$$X=(A*A^T)^{-1}A^TY$$

解出的X就是线性回归函数的系数矩阵。

边缘检测

Posted on 2016-02-21 | In Computer Vision

边缘检测在图像的检测中是经常会用到的。图片的边缘会包含大量的信息,因此在图像的分割、识别、分析中通常可以取边缘作为图像特征。边缘检测最经典的应用就是图像的锐化了,想必大家都用过。

为了进行边缘检测,我们通常会用到以下的一些算子,即一阶算子(梯度算子)和二阶算子(拉普拉斯算子)。

拉普拉斯算子(Laplacian Operator)

拉普拉斯算子其实就是二阶微分算子,具有各向同性。他也会增强噪声,因此在使用拉普拉斯算子之前需要进行平滑处理。他的主要用处其实是判断一个像素是处于图像的亮区还是暗区。

拉普拉斯算子定义式:
$$\begin{aligned}\nabla^2 f(x,y)&=\frac{\partial^2 f}{\partial x}+\frac{\partial^2 f}{\partial y}\\&=f(x+1,y)+f(x-1,y)+f(x,y+1)+f(x,y-1)-4f(x,y)\end{aligned}$$

常用模板公式:(1)$\begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}$(2)$\begin{bmatrix}1&1&1\\1&-8&1\\1&1&1\end{bmatrix}$

拉普拉斯模板有多个变种,上面是常用的两个,对于每一个像素,我们把他作为模板的中心元素,套用公式进行更新。需要注意的是,与梯度算子不同,他是各向同性的,因此不需要$G_x$,$G_y$两个模板,也就不能区分边缘的方向了。

还有一个常用的高斯-拉普拉斯算子(Gauss-Laplacian Operator),将平滑滤波器和拉普拉斯边缘检测算子结合起来的:

高斯-拉普拉斯算子模板:$\begin{bmatrix}-2&-4&-4&-4&-2\\-4&0&8&0&-4\\-4&8&24&8&-4\\-4&0&8&0&-4\\-2&-4&-4&-4&-2\end{bmatrix}$

另外,拉普拉斯算子有一个很实用的用途,就是对图像进行锐化,这时候只要将常用模板稍微改动一下即可得到锐化模板:

(1)$\begin{bmatrix}0&-1&0\\-1&5&-1\\0&-1&0\end{bmatrix}$(2)$\begin{bmatrix}-1&-1&-1\\-1&9&-1\\-1&-1&-1\end{bmatrix}$

测试

OK,最后就是用代码实现一下了,总体还是很简单的,而且类似于OpenCV之类的库都有实现的函数。不过我还是用pil来简单实现一下:

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
import Image,matplotlib,sys,ImageEnhance,ImageFont,ImageDraw
import numpy as np

def gauss(im_arr):
mat=np.array([[-2,-4,-4,-4,-2],[-4,0,8,0,-4],[-4,8,24,8,-4],[-4,0,8,0,4],[-2,-4,-4,-4,-2]])
height,width=im_arr.shape
ans=np.zeros((height-4,width-4))
for i in range(height-4):
for j in range(width-4):
ans[i,j]=np.sum(im_arr[i:i+5,j:j+5]*mat)
return ans

def laplacian(im_arr):
mat=np.array([[1,1,1],[1,-8,1],[1,1,1]])
height,width=im_arr.shape
ans=np.zeros((height-2,width-2))
for i in range(height-2):
for j in range(width-2):
ans[i,j]=np.sum(im_arr[i:i+3,j:j+3]*mat)
return ans

def sharpen(im_arr):
mat=np.array([[-1,-1,-1],[-1,9,-1],[-1,-1,-1]])
height,width=im_arr.shape
ans=np.zeros((height-2,width-2))
for i in range(height-2):
for j in range(width-2):
ans[i,j]=np.sum(im_arr[i:i+3,j:j+3]*mat)
return ans

def robert(im_arr):
height,width=im_arr.shape
ans=np.zeros((height-1,width-1))
for i in range(height-1):
for j in range(width-1):
ans[i,j]=np.abs(im_arr[i,j]-im_arr[i+1,j+1])+np.abs(im_arr[i+1,j]-im_arr[i,j+1])
return ans

def sobel(im_arr):
matx=np.array([[-1,0,1],[-2,0,2],[-1,0,1]])
maty=np.array([[1,2,1],[0,0,0],[-1,-2,-1]])
height,width=im_arr.shape
gx=np.zeros((height-2,width-2))
gy=np.zeros((height-2,width-2))
for i in range(height-2):
for j in range(width-2):
gx[i][j]=np.sum(matx*im_arr[i:i+3,j:j+3])
gy[i][j]=np.sum(maty*im_arr[i:i+3,j:j+3])
ans=np.abs(gx)+np.abs(gy)
return ans

def prewitt(im_arr):
matx=np.array([[1,0,-1],[1,0,-1],[1,0,-1]])
maty=np.array([[-1,-1,-1],[0,0,0],[1,1,1]])
height,width=im_arr.shape
gx=np.zeros((height-2,width-2))
gy=np.zeros((height-2,width-2))
for i in range(height-2):
for j in range(width-2):
gx[i][j]=np.sum(matx*im_arr[i:i+3,j:j+3])
gy[i][j]=np.sum(maty*im_arr[i:i+3,j:j+3])
ans=np.abs(gx)+np.abs(gy)
return ans

def isotropic(im_arr):
matx=np.array([[1,0,-1],[1.414,0,-1.414],[1,0,-1]])
maty=np.array([[-1,-1.414,-1],[0,0,0],[1,1.414,1]])
height,width=im_arr.shape
gx=np.zeros((height-2,width-2))
gy=np.zeros((height-2,width-2))
for i in range(height-2):
for j in range(width-2):
gx[i][j]=np.sum(matx*im_arr[i:i+3,j:j+3])
gy[i][j]=np.sum(maty*im_arr[i:i+3,j:j+3])
ans=np.abs(gx)+np.abs(gy)
return ans

def show(im_arr,name):
im=Image.fromarray(im_arr)
draw=ImageDraw.Draw(im)
font = ImageFont.truetype('/usr/share/fonts/truetype/liberation/LiberationMono-BoldItalic.ttf', 18)
draw.text((0,0),name,font=font,fill='red')
del draw
im=im.convert('L')
im.save(name+'.png')

if __name__=='__main__':
im=Image.open(sys.argv[1])
im=im.convert('L')
im.save('source.png')
show(laplacian(np.array(im,'int64')),'Laplace')
show(sobel(np.array(im,'int64')),'Sobel')
show(robert(np.array(im,'int')),'Robert')
show(prewitt(np.array(im,'int64')),'Prewitt')
show(isotropic(np.array(im,'int64')),'Isotropic')
show(sharpen(np.array(im,'int64')),'sharpen')
show(gauss(np.array(im,'int64')),'gauss-laplace')

test
source
Sobel
sharpen
Robert
Prewitt
Isotropic
Laplace
gauss-laplace

最后用他来处理一下lena.png,会的到上面的结果。

Julia集分形图案绘制

Posted on 2016-02-19 | In Others

从去年就开始窥东大的C++教学群,当时就被李骏扬老师讲的分形图案给吸引了,简直美赞了。他们的期末作业就是制作一个分形图案的视频,我们这种学校显然不会有这种东西。于是就想着能不能自己研究着画下,然而并不知道这种图案怎么画,度娘上找来的基本没用。搁置了一年,偶然间翻到了一篇论文,终于找到了画图的方法了,加上之前正好有用python绘图的工具,总算把这个东西搞通了一点。其实这个玩意的水还是非常深的,牵涉到了复分析,分形,甚至是混沌理论,据说从上古贝壳的图案,到如今麦田怪圈的图案,都和Julia集有关,说来也是玄乎。

Julia集

简单的讲,_Julia_集就是复平面上的一些点$z$,对于一个固定的复数$c$,$z$点在经过无限次$z\gets z^2+c$的迭代之后最终都收敛到一个固定的值上,那么复平面上所有这样的z点构成的集合就是_julia_集。(如果固定初始值而将$c$当做变量则生成的是mandelbrot集)

当然,这个迭代公式也有他的变种,比如多重julia集或者指数julia集等。

讲的不太清楚,具体的介绍可以参见_wiki pedia以及matrix67_的博客(对度娘表示无语),这里就不班门弄斧了。_wiki_里面主要介绍的是一些数学定义和推导以及他的一些典型图形,而_matrix67_写的则更加容易理解,他通过一步一步迭代过程的展现十分生动的描述了图像的产生过程。不过事实上,matrix67在之前的博客里虽然提供了绘图的代码,但是并没有介绍图像生成的算法。干看代码还是有点恶心的。而市面上又没怎么提及绘图的算法。事实上,他用的算法也非常的简单普遍,我们叫他“逃逸时间算法”。

逃逸时间算法(Escape Time Algorithm)

曾经纳闷了很长时间,上面那个简单的迭代式究竟是怎样生成那些纷繁复杂的分形图案的。最后终于在知网上找到了这个算法。

  1. 设定逃逸半径R,最高迭代次数N;
  2. 将初始的点z进行迭代,如果在N次迭代之内z的模超过了R,那么就认为z“逃逸出去”了,逃逸出去时的迭代次数n就是”逃逸时间“;
  3. 如果经过N次迭代,z的模仍然未到达R,那么就认为z是收敛集内的点;
  4. 将所有的逃逸点按照不同的逃逸时间进行染色就得到了美丽的“逃逸时间图”。

其实逃逸时间图显示的并不是真正意义上的_julia_集,而是不属于_julia_集合的点。

当然,还有一种常用的_julia_集绘图算法–外部距离估计算法,这里不做过多介绍。

静态实现

用python+matplotlib模拟逃逸时间算法进行简单绘图:

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
import numpy as np
import matplotlib.pyplot as plt
import time

c_real=-1.621
c_imag=-2

def func(a):
c=np.complex(c_real,c_imag)
return np.exp(a*a*a)+c+a*a-a

R=50
density=400
N=100

x=[]
y=[]
c=[]

for i in xrange(density):
for j in xrange(density):
real=(i-density/2.0)/density*2.0
imag=(j-density/2.0)/density*2.0
now=np.complex(real,imag)
for n in xrange(N):
now=func(now)
if np.abs(now)>R:
x.append(real)
y.append(imag)
c.append(n+1)
break

plt.scatter(x,y,s=2,c=c,edgecolors='face')
plt.xticks([])
plt.yticks([])
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.savefig(str(const_real)+'+'+str(const_imag)+'i.png')

参数c由c_real和c_imag确定,迭代公式由func函数确定,描点时注意去掉边界以及标尺。

通过这个模板可以非常轻松的绘出下面这几张图:

c=-0.74543+0.11301i

c=-0.8+0.156i

c=-0.4+0.6i

c=-0.621+0i,func=exp(z*z*z)+c

c=0.285+0.01i

动态实现

既然我们可以直接绘出迭代N次之后的图像,那么我们当然也可以画出每一次的图像组合成动图,这样看起来效果也更明显。

图像生成:

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
50
51
52
53
54
55
56
57
58
import numpy as np
import matplotlib.pyplot as plt
import time

def func(a):
c=np.complex(-0.74543,0.11301)
return a*a+c

R=10
density=400
N=100

time1=time.time()

mat=[]
for i in xrange(density):
for j in xrange(density):
real=(i-density/2.0)/density*2.0
imag=(j-density/2.0)/density*2.0
now=np.complex(real,imag)
for n in xrange(N):
now=func(now)
if np.abs(now)>R:
mat.append((real,imag,n+1))
break

time2=time.time()

print 'calculating costs '+str(time2-time1)+' s.'

mat.sort(lambda x,y:cmp(x[2],y[2]))

time3=time.time()
print 'sorting costs '+str(time3-time2)+' s.'

real=[]
imag=[]
color=[]

cnt=1
for i in xrange(len(mat)-1):
real.append(mat[i][0])
imag.append(mat[i][1])
color.append(mat[i][2])
if i==len(mat)-2 or (not mat[i][2]==mat[i+1][2]):
plt.scatter(real,imag,s=2,c=color,alpha=1,edgecolors='face')
plt.xticks([])
plt.yticks([])
plt.xlim(-1,1)
plt.ylim(-1,1)
plt.savefig('%02d' % cnt)
print 'pic'+('%02d' % cnt)+' OK.'
cnt+=1

time4=time.time()
print 'generating costs '+str(time4-time3)+' s.'

print 'total costs '+str(time4-time3)+' s.'

生成图像的过程比较缓慢,相比而言计算所耗费的时间就不算太多了。

gif组合:

1
$convert *.png -resize 320x240 out.gif

我们用ImageMagick工具来生成gif,这里如果用-layer来压缩会导致图像失真,因此用改变图像大小的方式来适当压缩文件。

最后就可以生成动图:

Ubuntu下Opencv安装使用简述

Posted on 2016-02-11 | In Computer Vision

简述

Opencv就不解释了,是个很有名的图形库。不仅在进行软件开发的过程中需要用到,而且他也是很多开源软件的运行依赖,所以安装一个Opencv就很有必要了,即使自己本身并不想学习使用。

安装

以下主要是从百度上找到的可用方法:

安装运行依赖

1
2
3
4
5
$ sudo apt-get install libqt4-dev libopencv-dev build-essential cmake git libgtk2.0-dev pkg-config\
python-dev python-numpy libdc1394-22 libdc1394-22-dev libjpeg-dev libpng12-dev libtiff5-dev \
libjasper-dev libavcodec-dev libavformat-dev libswscale-dev libxine2-dev libgstreamer0.10-dev\
libgstreamer-plugins-base0.10-dev libv4l-dev libtbb-dev libfaac-dev libmp3lame-dev libopencore-amrnb-dev \
libopencore-amrwb-dev libtheora-dev libvorbis-dev libxvidcore-dev x264 v4l-utils unzip

下载源代码

在 官网 下载适合的版本就好,我这下的是3.1.0。

下载完成解压就好。

编译

编译还是有点麻烦的,现在都是用cmake结合make来编译,头一次用这个有点头大。

  1. 在文件目录下新建一个叫build/的文件夹,这个文件夹的目的主要是用来存放编译生成的临时文件。当然起别的名字也可以。
  2. 进入build/文件夹,输入以下命令
    1
    cmake -D CMAKE_BUILD_TYPE=RELEASE -D CMAKE_INSTALL_PREFIX=/usr/local  -D BUILD_PYTHON_SUPPORT=ON       -D WITH_XINE=ON       -D WITH_OPENGL=ON       -D WITH_TBB=ON       -D BUILD_EXAMPLES=ON       -D BUILD_NEW_PYTHON_SUPPORT=ON       -D WITH_V4L=ON  ..
    其实就是cmake命令后加些配置参数,最后是CMakeLists.txt配置文件的位置,这个位置当然就是../
  3. 编译安装,输入make -j $(nproc)(这是多进程make的命令,$(nproc)就是进程数,当然也可以直接指定),这会花很长时间,最后再$ sudo make install ,对文件进行安装。
  4. 最后还要配置一些路径,输入以下命令
    1
    /bin/bash -c 'echo "/usr/local/lib" > /etc/ld.so.conf.d/opencv.conf'
    然后再ldconfig即可。(注意给权限)
  5. 最后可以用如下命令来判断是否安装成功
    1
    2
    pkg-config --modversion opencv
    pkg-config --cflags opencv

测试

OpenCV在codeblock下可以通过加链接库的形式编译运行。但是如果在命令行下,就得手写cmake了。

为了方便测试,我们新建一个test文件夹,在这下面写一个测试程序。

首先新建如下文件,保存为test.cpp

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
#include<opencv2/highgui.hpp>
#include<opencv2/imgproc.hpp>
using namespace cv;
int main(int argc ,char** argv){
if(argc!=2){
printf("No image data\n");
return -1;
}
char *imageName=argv[1];
Mat image;
image=imread(imageName,1);
if(!image.data){
printf("No iamge data\n");
return -1;
}
namedWindow(imageName,CV_WINDOW_AUTOSIZE);
imshow(imageName,image);
waitKey(0);
return 0;
}

然后随便把一个测试用图片复制到test文件夹下,我用的是他自带的最经典的lena.jpg。

接着编写cmake配置文件,将下面的文件保存为CMakeLists.txt

1
2
3
4
project(test)
add_executable(test test.cpp)
find_package(OpenCV REQUIRED)
target_link_libraries(test ${OpenCV_LIBS})

里面用到的各种文件的意义也很清楚了,以后照着改就行。

最后在test下新建build文件夹,进入后输入cmake ../ 即可完成cmake编译,然后再输入make即可生成可执行文件。

找到test文件,然而在命令行下输入./test ../lena.jpg 即可运行程序。

libsvm库的使用

Posted on 2016-02-06 | In Machine Learning

看了下svm(支持向量机)的实现原理,感觉基础的部分还是不难懂的,但是如果要自己动手实现的话还是有很大难度的,况且自己写的效果肯定不太好。于是就在网上找了一个大牛写的svm库,实现了多种分类方式,而且涵盖了几乎所有常见语言的接口,用起来方便而且效果也很好。

概述

LIBSVM是台湾大学林智仁(Lin Chih-Jen)教授等开发设计的,综合使用了包括线性函数,多项式函数,径向基函数,sigmoid函数等在内的不同分类方式,而且支持包括C/C++,python,java,matlab,Octave,R,C#,Perl,Ruby,Weka,Node.js,Scilab,Lisp,,haskell,Cuba,.Net,PHP等等一系列语言(听上去就很强大)。

具体的介绍可以参见林智仁教授的HomePage,这里面有他个人的介绍,以及Libsvm的各种文档以及最新版本。

安装

最好的办法就是去林教授的主页上下了,可以见到最完整的文件以及文档,而且还附带基础教程的测试样例。不过调用起来不太方便,还得手动将源文件配置到正确的地方才能随时使用。

我为了图省事,就直接用apt下载了(没想到ubuntu竟然收录了这个,可见这玩意的强大)。

1
$sudo apt-get install python-libsvm

这样我们就能像平常导入包一样在python中调用了。

当然我们需要一些文档,试着用man来查看发现并没有,于是locate了一下发现他的文档在这里:

1
2
/usr/share/doc/libsvm3/README.gz
/usr/share/doc/python-libsvm/README.gz

OK,需要用到的时候进去看一下就好了。

基本用法

教程中对libsvm的用法有两个档次,即 “high-level ”和”low-level” 。这里的 high-level 并不是指高端用法,low-level 也并不是指低端用法(话说我一开始就是这么理解的0.0)。其实这里的 high-level 是指封装程度高,也就是细节隐藏的更好,用户使用更方便;同样,low-level 是指所用的函数更加底层,更加体现细节,但是用起来难度就更高了。当然,我会选择high-level来进行应用。

具体函数的用法都在文档里,就不一一记录了,我这里姑且就对他进行一个简要的应用。下面就用libsvm来代替之前在 正方系统验证码识别 项目中的那个logistic_sgd.py文件,并顺便查看下svm算法的效果。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
import svmutil,cPickle,gzip
import numpy as np

f = gzip.open('vericode.pkl.gz')
train_set,valid_set,test_set = cPickle.load(f)
f.close()

x = train_set[0].tolist()
x.extend(valid_set[0].tolist())

y = train_set[1].tolist()
y.extend(valid_set[1].tolist())

#以上是为了控制输入格式,y是用list储存的分类标签;x是用list储存的数据,每一个数据是一个list,保存该数据的特征值

prob = svmutil.svm_problem(y,x)#用数据生成表示训练集的对象

param = svmutil.svm_parameter('-t 0')#表示训练参数,这里选择的是线性分类

#m = svm_load_model('heart_scale.model')可以导入之前得到的最优模型
m = svmutil.svm_train(prob,param)#开始训练,并生成训练后的模型
svmutil.svm_predict(y,x, m)#检测训练集的符合程度
p_label, p_acc, p_val = svmutil.svm_predict(test_set[1].tolist(), test_set[0].tolist(), m)#检测测试集的泛化误差,当然也可以用来进行预测,p_label就是用来保存预测结果的
#svmutil.svm_save_model('heart_scale.model', m)#可以将最优模型保存下来

本地运行结果:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
*
optimization finished, #iter = 37
nu = 0.003704
obj = -0.251887, rho = 0.269514
nSV = 10, nBSV = 0
*
optimization finished, #iter = 54
nu = 0.010331
obj = -0.687010, rho = -0.920502
nSV = 16, nBSV = 0
......
*
optimization finished, #iter = 96
nu = 0.004130
obj = -0.280826, rho = 0.255170
nSV = 21, nBSV = 0
*.*
optimization finished, #iter = 178
nu = 0.009597
obj = -0.767762, rho = -0.705452
nSV = 29, nBSV = 0
Total nSV = 1409
Accuracy = 99.009% (2198/2220) (classification)
Accuracy = 91.3669% (508/556) (classification)

前面一大段的内容是在svm_train训练的时候自动输出的,大概反映了训练的进程(具体意义暂时没搞懂)。

最后两个分别是对训练集的拟合程度,以及对测试集的泛化拟合程度。我们可以看见他的拟合效果非常的好,几乎达到了完美。后来的泛化拟合效果稍微有点问题,这大概是出现了过拟合问题,而这个问题通常是可以通过加大训练集的量来解决的。总的来说,训练的效率还是非常高的,用起来也非常的方便。

自制正方软件系统验证码的识别程序(4/4)

Posted on 2016-02-05 | In Computer Vision

效果

最后总的效果还是不错的,从测试集上面看,单个字符的识别准确率平均能达到90%左右,最好的时候达到过93%。我想这也差不多达到了优化的极限了,毕竟在图片样子不标准,位置更不标准。而且仅仅在截取子图这一步上就会与实际情况有些出入,况且一些细微的差别就算是人脑也容易出错,比如1和7、i和j等等。加上我使用整个数据集的量不算大,满打满算也就500张图,所以能进行33个字符的识别已经挺不错的了。官网上60000张图的数据集的学习也差不多达到了93%而已。

总结

代码就是这样了,不算多,也就小几百行。但是很多细节的把握以及参数的调整却让我头大了很久。

开始的时候也只是照着mnist的样子直接二值化后(仅仅调用了blur1函数)交给学习机来学习,效果可想而知。。。正确率只有可怜的10%+,让我一度怀疑这个算法的可行性。然而冷静下来想了想是不是因为线条过细导致信息量不够大。于是就像办法加粗线条,最后选择了用高斯模糊来进行边缘处理再二值化。果然,这么一搞就把正确率调到了50%左右。但是这个识别率也实在是太低了,想了半天又没啥办法,试着调节学习算法中的参数效果也并不理想。最后盯着8的两个圈圈看的时候突然想起来是不是可以通过数圈圈的个数来区分一些数字!于是写了一个种子填充,果然将正确率提高到了70%。有了这个经验,我想都没想就再加进去了count_fill和count_border函数计算两个特征值,最后将正确率提高到了90%。

除了算法设计部分,由于是第一次使用PIL,对图片的处理让我蛋疼了好久,比如二值的图像和灰度的图像之间性质的差别问题,图像矩阵中int到bool的转化问题,内存中的图像对象的显示错误问题等等,都特别让人伤脑筋。(不过多谢度娘)

还有,这次也是我第一次用python写这种小项目。深刻觉得python这种脚本语言还真是方便,比如ipython的即时反馈,数组列表字典元组的使用,文件和文件夹的IO,函数的调用和使用(竟然可以在函数的内部定义另一个函数,一个函数就能做到C中一个类才能做到的事情)等等。不过调试起来还是有点麻烦的,万一程序有点小错,他报的错经常会跳到好多地方去,让人摸不着头脑。

虽然这种代码的意义不大(毕竟只适用于这个特定的验证码),但是做为自己写的第一个能用的机器学习程序,还是挺有里程碑意义的。自己mark下。

附github上的完整项目:Zhengfang_CaptchaRecognition

自制正方软件系统验证码的识别程序(3/4)

Posted on 2016-02-05 | In Computer Vision

util.py

这个文件里主要提供了5个函数,提供给package.py使用,特别是对特征值的计算。

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
import Image,os,ImageFilter
import numpy as np
table1=[]
table2=[]
threshold1=18
for i in range(256):
if i<threshold1:
table1.append(0)
else:
table1.append(1)
threshold2=240
for i in range(256):
if i<threshold2:
table2.append(0)
else:
table2.append(1)

def blur1(im):
im=im.convert('L')
im=im.point(table1,'1')
return im
def blur2(im):
im=im.convert('L')
im=im.filter(ImageFilter.GaussianBlur(1.5))
im=im.point(table2,'1')
return im

def seed_fill(im):
arr=np.array(im.convert('L')).astype('bool').astype('int')
height,width=arr.shape
arr=arr.tolist()

output=[]
output.append(np.ones(width+2,).astype('int').tolist())
for i in range(height):
tmp=[]
tmp.append(1)
for j in range(width):
tmp.append(arr[i][j])
tmp.append(1)
output.append(tmp)
output.append(np.ones(width+2,).astype('int').tolist())
arr=output
height+=2
width+=2
def dfs(i,j):
if(arr[i][j]==0):
return 0
arr[i][j]=0
if(i>0):
dfs(i-1,j)
if(i<height-1):
dfs(i+1,j)
if(j>0):
dfs(i,j-1)
if(j<width-1):
dfs(i,j+1)
return 1
cnt=0
for i in range(height):
for j in range(width):
cnt+=dfs(i,j)
return cnt

def count_border(im):
arr=np.array(im.convert('L')).astype('bool').astype('int')
height,width=arr.shape
cnt=0
for i in range(height):
for j in range(width):
if arr[i][j]==1:
continue
if i==0 or i==height-1 or j==0 or j==width-1 :
cnt+=1
continue
if arr[i-1][j]==1 or arr[i+1][j]==1 or arr[i][j-1]==1 or arr[i][j+1]==1:
cnt+=1
return cnt

def count_fill(im):
arr=np.array(im.convert('L')).astype('bool').astype('int')
height,width=arr.shape
cnt=0
for i in range(height):
for j in range(width):
if arr[i][j]==0:
cnt+=1
return cnt
  • blur1函数实现的是对图片先进行灰度化,再进行二值化,结合table1数组的使用,将亮度大于某一较低阈值的统统设为纯白,否则为纯黑。这样就只保留了字符的骨架,生成的图片十分便于特征值的测量与计算,而且极大消除了噪声的影响。
  • blur2函数实现的是对图片的模糊化再二值化,为的是将图片变“粗”一点,便于后续的学习(这也是后来偶然发现的提高学习准确性的方法)。注意模糊参数的调节。
  • seed_fill函数实现种子填充,即返回空白联通块的个数,比如8返回3,0返回2等等。当然,首先得在外围加一圈虚拟的白圈防止边界的阻隔。
  • count_border和count_fill函数实现的是边界黑点的计数,以及黑点数目的计数。这只是我开脑洞想到的特征,好像也有点用的。

logistic_sgd.py

这是官网上的学习的核心算法,我只是稍微修改了下参数还有predict函数,有些实现细节还没有彻底搞懂(话说theano还真不怎么好理解)。

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
"""
This tutorial introduces logistic regression using Theano and stochastic
gradient descent.

Logistic regression is a probabilistic, linear classifier. It is parametrized
by a weight matrix :math:`W` and a bias vector :math:`b`. Classification is
done by projecting data points onto a set of hyperplanes, the distance to
which is used to determine a class membership probability.

Mathematically, this can be written as:

.. math::
P(Y=i|x, W,b) &= softmax_i(W x + b) \\
&= \frac {e^{W_i x + b_i}} {\sum_j e^{W_j x + b_j}}

The output of the model or prediction is then done by taking the argmax of
the vector whose i'th element is P(Y=i|x).

.. math::

y_{pred} = argmax_i P(Y=i|x,W,b)

This tutorial presents a stochastic gradient descent optimization method
suitable for large datasets.

References:

- textbooks: "Pattern Recognition and Machine Learning" -
Christopher M. Bishop, section 4.3.2

"""
__docformat__ = 'restructedtext en'

import cPickle
import matplotlib.pyplot as plt
import gzip
import os
import sys
import timeit
import numpy
import numpy as np

import theano
import theano.tensor as T

class LogisticRegression(object):
"""Multi-class Logistic Regression Class

The logistic regression is fully described by a weight matrix :math:`W`
and bias vector :math:`b`. Classification is done by projecting data
points onto a set of hyperplanes, the distance to which is used to
determine a class membership probability.
"""

def __init__(self, input, n_in, n_out):
""" Initialize the parameters of the logistic regression

:type input: theano.tensor.TensorType
:param input: symbolic variable that describes the input of the
architecture (one minibatch)

:type n_in: int
:param n_in: number of input units, the dimension of the space in
which the datapoints lie

:type n_out: int
:param n_out: number of output units, the dimension of the space in
which the labels lie

"""
# start-snippet-1
# initialize with 0 the weights W as a matrix of shape (n_in, n_out)
self.W = theano.shared(
value=numpy.zeros(
(n_in, n_out),
dtype=theano.config.floatX
),
name='W',
borrow=True
)
# initialize the biases b as a vector of n_out 0s
self.b = theano.shared(
value=numpy.zeros(
(n_out,),
dtype=theano.config.floatX
),
name='b',
borrow=True
)

# symbolic expression for computing the matrix of class-membership
# probabilities
# Where:
# W is a matrix where column-k represent the separation hyperplane for
# class-k
# x is a matrix where row-j represents input training sample-j
# b is a vector where element-k represent the free parameter of
# hyperplane-k
self.p_y_given_x = T.nnet.softmax(T.dot(input, self.W) + self.b)

# symbolic description of how to compute prediction as class whose
# probability is maximal
self.y_pred = T.argmax(self.p_y_given_x, axis=1)
# end-snippet-1

# parameters of the model
self.params = [self.W, self.b]

# keep track of model input
self.input = input

def negative_log_likelihood(self, y):
"""Return the mean of the negative log-likelihood of the prediction
of this model under a given target distribution.

.. math::

\frac{1}{|\mathcal{D}|} \mathcal{L} (\theta=\{W,b\}, \mathcal{D}) =
\frac{1}{|\mathcal{D}|} \sum_{i=0}^{|\mathcal{D}|}
\log(P(Y=y^{(i)}|x^{(i)}, W,b)) \\
\ell (\theta=\{W,b\}, \mathcal{D})

:type y: theano.tensor.TensorType
:param y: corresponds to a vector that gives for each example the
correct label

Note: we use the mean instead of the sum so that
the learning rate is less dependent on the batch size
"""
# start-snippet-2
# y.shape[0] is (symbolically) the number of rows in y, i.e.,
# number of examples (call it n) in the minibatch
# T.arange(y.shape[0]) is a symbolic vector which will contain
# [0,1,2,... n-1] T.log(self.p_y_given_x) is a matrix of
# Log-Probabilities (call it LP) with one row per example and
# one column per class LP[T.arange(y.shape[0]),y] is a vector
# v containing [LP[0,y[0]], LP[1,y[1]], LP[2,y[2]], ...,
# LP[n-1,y[n-1]]] and T.mean(LP[T.arange(y.shape[0]),y]) is
# the mean (across minibatch examples) of the elements in v,
# i.e., the mean log-likelihood across the minibatch.
return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]), y])
# end-snippet-2

def errors(self, y):
"""Return a float representing the number of errors in the minibatch
over the total number of examples of the minibatch ; zero one
loss over the size of the minibatch

:type y: theano.tensor.TensorType
:param y: corresponds to a vector that gives for each example the
correct label
"""

# check if y has same dimension of y_pred
if y.ndim != self.y_pred.ndim:
raise TypeError(
'y should have the same shape as self.y_pred',
('y', y.type, 'y_pred', self.y_pred.type)
)
# check if y is of the correct datatype
if y.dtype.startswith('int'):
# the T.neq operator returns a vector of 0s and 1s, where 1
# represents a mistake in prediction
return T.mean(T.neq(self.y_pred, y))
else:
raise NotImplementedError()

def load_data(dataset):
''' Loads the dataset

:type dataset: string
:param dataset: the path to the dataset (here MNIST)
'''

#############
# LOAD DATA #
#############
''''
# Download the MNIST dataset if it is not present
data_dir, data_file = os.path.split(dataset)
if data_dir == "" and not os.path.isfile(dataset):
# Check if dataset is in the data directory.
new_path = os.path.join(
os.path.split(__file__)[0],
"..",
"data",
dataset
)
if os.path.isfile(new_path) or data_file == 'mnist.pkl.gz':
dataset = new_path

if (not os.path.isfile(dataset)) and data_file == 'mnist.pkl.gz':
import urllib
origin = (
'http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz'
)
print 'Downloading data from %s' % origin
urllib.urlretrieve(origin, dataset)
'''
print '... loading data'

# Load the dataset
f = gzip.open(dataset, 'rb')
train_set, valid_set, test_set = cPickle.load(f)
f.close()
#train_set, valid_set, test_set format: tuple(input, target)
#input is an numpy.ndarray of 2 dimensions (a matrix)
#witch row's correspond to an example. target is a
#numpy.ndarray of 1 dimensions (vector)) that have the same length as
#the number of rows in the input. It should give the target
#target to the example with the same index in the input.

def shared_dataset(data_xy, borrow=True):
""" Function that loads the dataset into shared variables

The reason we store our dataset in shared variables is to allow
Theano to copy it into the GPU memory (when code is run on GPU).
Since copying data into the GPU is slow, copying a minibatch everytime
is needed (the default behaviour if the data is not in a shared
variable) would lead to a large decrease in performance.
"""
data_x, data_y = data_xy
shared_x = theano.shared(numpy.asarray(data_x,
dtype=theano.config.floatX),
borrow=borrow)
shared_y = theano.shared(numpy.asarray(data_y,
dtype=theano.config.floatX),
borrow=borrow)
# When storing data on the GPU it has to be stored as floats
# therefore we will store the labels as ``floatX`` as well
# (``shared_y`` does exactly that). But during our computations
# we need them as ints (we use labels as index, and if they are
# floats it doesn't make sense) therefore instead of returning
# ``shared_y`` we will have to cast it to int. This little hack
# lets ous get around this issue
return shared_x, T.cast(shared_y, 'int32')

test_set_x, test_set_y = shared_dataset(test_set)
valid_set_x, valid_set_y = shared_dataset(valid_set)
train_set_x, train_set_y = shared_dataset(train_set)

rval = [(train_set_x, train_set_y), (valid_set_x, valid_set_y),
(test_set_x, test_set_y)]
return rval

def sgd_optimization_mnist(learning_rate=0.1, n_epochs=10000,
dataset='vericode.pkl.gz',
batch_size=100):
"""
Demonstrate stochastic gradient descent optimization of a log-linear
model

This is demonstrated on MNIST.

:type learning_rate: float
:param learning_rate: learning rate used (factor for the stochastic
gradient)

:type n_epochs: int
:param n_epochs: maximal number of epochs to run the optimizer

:type dataset: string
:param dataset: the path of the MNIST dataset file from
http://www.iro.umontreal.ca/~lisa/deep/data/mnist/mnist.pkl.gz

"""
datasets = load_data(dataset)

train_set_x, train_set_y = datasets[0]
valid_set_x, valid_set_y = datasets[1]
test_set_x, test_set_y = datasets[2]

# compute number of minibatches for training, validation and testing
n_train_batches = train_set_x.get_value(borrow=True).shape[0] / batch_size
n_valid_batches = valid_set_x.get_value(borrow=True).shape[0] / batch_size
n_test_batches = test_set_x.get_value(borrow=True).shape[0] / batch_size

######################
# BUILD ACTUAL MODEL #
######################
print '... building the model'

# allocate symbolic variables for the data
index = T.lscalar() # index to a [mini]batch

# generate symbolic variables for input (x and y represent a
# minibatch)
x = T.matrix('x') # data, presented as rasterized images
y = T.ivector('y') # labels, presented as 1D vector of [int] labels

# construct the logistic regression class
classifier = LogisticRegression(input=x, n_in=12 * 18+3, n_out=36)

# the cost we minimize during training is the negative log likelihood of
# the model in symbolic format
cost = classifier.negative_log_likelihood(y)

# compiling a Theano function that computes the mistakes that are made by
# the model on a minibatch
test_model = theano.function(
inputs=[index],
outputs=classifier.errors(y),
givens={
x: test_set_x[index * batch_size: (index + 1) * batch_size],
y: test_set_y[index * batch_size: (index + 1) * batch_size]
}
)

validate_model = theano.function(
inputs=[index],
outputs=classifier.errors(y),
givens={
x: valid_set_x[index * batch_size: (index + 1) * batch_size],
y: valid_set_y[index * batch_size: (index + 1) * batch_size]
}
)

# compute the gradient of cost with respect to theta = (W,b)
g_W = T.grad(cost=cost, wrt=classifier.W)
g_b = T.grad(cost=cost, wrt=classifier.b)

# start-snippet-3
# specify how to update the parameters of the model as a list of
# (variable, update expression) pairs.
updates = [(classifier.W, classifier.W - learning_rate * g_W),
(classifier.b, classifier.b - learning_rate * g_b)]

# compiling a Theano function `train_model` that returns the cost, but in
# the same time updates the parameter of the model based on the rules
# defined in `updates`
train_model = theano.function(
inputs=[index],
outputs=cost,
updates=updates,
givens={
x: train_set_x[index * batch_size: (index + 1) * batch_size],
y: train_set_y[index * batch_size: (index + 1) * batch_size]
}
)
# end-snippet-3

###############
# TRAIN MODEL #
###############
print '... training the model'
# early-stopping parameters
patience = 5000 # look as this many examples regardless
patience_increase = 2 # wait this much longer when a new best is
# found
improvement_threshold = 0.995 # a relative improvement of this much is
# considered significant
validation_frequency = min(n_train_batches, patience / 2)
# go through this many
# minibatche before checking the network
# on the validation set; in this case we
# check every epoch

best_validation_loss = numpy.inf
test_score = 0.
start_time = timeit.default_timer()

done_looping = False
epoch = 0
while (epoch < n_epochs) and (not done_looping):
epoch = epoch + 1
for minibatch_index in xrange(n_train_batches):
minibatch_avg_cost = train_model(minibatch_index)
# iteration number
iter = (epoch - 1) * n_train_batches + minibatch_index

if (iter + 1) % validation_frequency == 0:
# compute zero-one loss on validation set
validation_losses = [validate_model(i)
for i in xrange(n_valid_batches)]
this_validation_loss = numpy.mean(validation_losses)

print(
'epoch %i, minibatch %i/%i, validation error %f %%' %
(
epoch,
minibatch_index + 1,
n_train_batches,
this_validation_loss * 100.
)
)

# if we got the best validation score until now
if this_validation_loss < best_validation_loss:
#improve patience if loss improvement is good enough
if this_validation_loss < best_validation_loss * \
improvement_threshold:
patience = max(patience, iter * patience_increase)

best_validation_loss = this_validation_loss
# test it on the test set

test_losses = [test_model(i)
for i in xrange(n_test_batches)]
test_score = numpy.mean(test_losses)

print(
(
' epoch %i, minibatch %i/%i, test error of'
' best model %f %%'
) %
(
epoch,
minibatch_index + 1,
n_train_batches,
test_score * 100.
)
)

# save the best model
with open('best_model.pkl', 'w') as f:
cPickle.dump(classifier, f)

if patience <= iter:
done_looping = True
break

end_time = timeit.default_timer()
print(
(
'Optimization complete with best validation score of %f %%,'
'with test performance %f %%'
)
% (best_validation_loss * 100., test_score * 100.)
)
print 'The code run for %d epochs, with %f epochs/sec' % (
epoch, 1\. * epoch / (end_time - start_time))
print >> sys.stderr, ('The code for file ' +
os.path.split(__file__)[1] +
' ran for %.1fs' % ((end_time - start_time)))

def predict():
"""
An example of how to load a trained model and use it
to predict labels.
"""

# load the saved model
classifier = cPickle.load(open('best_model.pkl'))

# compile a predictor function
predict_model = theano.function(
inputs=[classifier.input],
outputs=classifier.y_pred)

# We can test it on some examples from test test
dataset='vericode.pkl.gz'
datasets = cPickle.load(gzip.open(dataset))
test_set_x, test_set_y = datasets[2]
predicted_values = predict_model(test_set_x[:10])
test_set_y=test_set_y[:10].astype('int')
dict='012345678abcdefghijklmnpqrstuvwxy'
pred=[]
ans=[]
for i in range(len(predicted_values)):
pred.append(dict[predicted_values[i]])
ans.append(dict[test_set_y[i]])
print ("Predicted values for the first 10 examples in test set:")
print pred
print ("Correct values for the first 10 examples in test set:")
print ans

if __name__ == '__main__':
sgd_optimization_mnist()
predict()

主要需要关注的就是 sgd_optimization_mnist 函数的传参,比如学习因子$\alpha$还有分组的大小。迭代的次数反倒不用管他,调到最高也成。

train.py

训练接口,很简单0.0,只是为了方便操作而已。

1
2
3
4
5
import os
os.system('python split.py')
os.system('python util.py')
os.system('python package.py')
os.system('python logistic_sgd.py')

运行之前所有的程序,然后直接训练。split.py会生成名叫number/的文件加,package.py会生成名叫vericode.pkl.gz的数据集,而logistic_sgd.py会直接进行训练,并把最好的结果输出到** best_model.pkl** 中加以保存,并在命令行中打印类似下面的信息:

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
... spliting
... packaging
... loading data
... building the model
... training the model
epoch 1, minibatch 16/16, validation error 71.200000 %
epoch 1, minibatch 16/16, test error of best model 70.800000 %
epoch 2, minibatch 16/16, validation error 56.800000 %
epoch 2, minibatch 16/16, test error of best model 57.600000 %
epoch 3, minibatch 16/16, validation error 49.600000 %
epoch 3, minibatch 16/16, test error of best model 50.200000 %
epoch 4, minibatch 16/16, validation error 44.200000 %
epoch 4, minibatch 16/16, test error of best model 43.000000 %
epoch 5, minibatch 16/16, validation error 38.800000 %
epoch 5, minibatch 16/16, test error of best model 40.000000 %
epoch 6, minibatch 16/16, validation error 36.400000 %
......
epoch 678, minibatch 16/16, validation error 10.000000 %
epoch 679, minibatch 16/16, validation error 10.000000 %
epoch 680, minibatch 16/16, validation error 10.200000 %
epoch 681, minibatch 16/16, validation error 10.200000 %
Optimization complete with best validation score of 9.600000 %,with test performance 11.400000 %
The code run for 682 epochs, with 60.733475 epochs/sec
Predicted values for the first 10 examples in test set:
['l', 'j', 'r', '0', 'q', '3', '6', '3', '7', 'p']
Correct values for the first 10 examples in test set:
['l', 'j', 'r', '0', 'q', '3', '6', '3', '7', 'p']

细致的显示了训练集的训练结果。

check.py

这个文件封装了调用训练成果进行识别的接口。由于懒得优化之前的代码结构,这里没有实现代码重用而是相当于重新写了一边之前的图片处理的过程。

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
import cPickle,os,Image,ImageFilter,theano,sys,time
import numpy as np
from logistic_sgd import LogisticRegression
from util import *

def recognize(pic='test'):
im=Image.open(pic)

#split
image=[]
image.append(im.crop((5,2,17,20)))
image.append(im.crop((17,2,29,20)))
image.append(im.crop((29,2,41,20)))
image.append(im.crop((41,2,53,20)))

#package
arr=[]
for img in image:
img=blur1(img)
cnt1=seed_fill(img)
cnt2=count_fill(img)
cnt2=(cnt2-50)/50
cnt3=count_border(img)
cnt3=(cnt3-50)/50
img=blur2(img)
array=np.array(img)
array=array.reshape(12*18).astype('int').astype('bool').astype('float32')
array=array.tolist()
array.append(cnt1)
array.append(cnt2)
array.append(cnt3)
array=np.array(array).astype('float32')
arr.append(array)

arr=np.array(arr)

#predict
classifier = cPickle.load(open('best_model.pkl'))
predict_model = theano.function(inputs=[classifier.input],outputs=classifier.y_pred)
predicted_values = predict_model(arr)
dict='012345678abcdefghijklmnpqrstuvwxy'
pred=[]
for i in range(len(predicted_values)):
pred.append(dict[predicted_values[i]])
return pred[0]+pred[1]+pred[2]+pred[3]

def recur_recognize(path):
if not path[-1]=='/':
path.join('/')
files=os.listdir(path)
cost=time.time()
for i in files:
#print '.',
#sys.stdout.flush()
old_name=path+str(i)
name=recognize(old_name)
new_name=path+str(name)
os.rename(old_name,new_name)

return time.time()-cost,len(files)

if __name__ == '__main__':
if len(sys.argv) <2:
print 'Please enter the file'
else:
if os.path.isdir(sys.argv[1]):
cost,size= recur_recognize(sys.argv[1])
print '\nRecognized '+str(size)+' pictures , in '+str(cost)+' s\n'
else:
print recognize(sys.argv[1])

这里实现了对单个文件的识别以及文件夹下所有文件的批量识别。可以如下使用:

$python check.py test.png 以及$python check.py png/

自制正方软件系统验证码的识别程序(2/4)

Posted on 2016-02-05 | In Computer Vision

文件组成

为了实现训练以及识别的过程,我总共设计了6个文件,作用如下:

文件 作用
split.py 用于将验证码中四个小字符分割出来,并分类保存。
util.py 用于保存一些常用的函数
logistic_sgd.py 这是官网上的样例代码,实现了softmax的分类算法,当然还要进行下修改
package.py 这个用来将图像数据进行处理并打包压缩成为方便使用的数据集
train.py 这是开始训练的接口
check.py 这是利用训练结果进行识别的接口

还有两个文件夹:

文件夹 作用
recognized/ 用于保存可作为训练集的图片,图片以内容的实际值命名
number/ 用于保存分离子图后的小图片,以类似0/ 1/ a/….的文件夹的形式进行分类

运行依赖

这里我使用的是python2.7版,需要以下必要的运行库的支持:

theano:

安装: $sudo pip install theano

目的:这是学习算法的运行依赖

numpy:

安装:$sudo apt-get install python-numpy

目的:提供了很多必要的数学运算方法

PIL:

安装:$sudo apt-get install python-pil

目的:进行图片的处理

split.py

由于牵涉到文件写入,而且与其他代码不存在依赖关系,所以把他独立出来方便修改。

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
import os
import numpy as np
import Image
import shutil

print '... spliting'

source_dir='recognized/'
dest_dir='number/'

if os.path.exists(dest_dir):
shutil.rmtree(dest_dir)

os.mkdir(dest_dir)
list=os.listdir(source_dir)
directory='123456780abcdefghijklmnpqrstuvwxy'
cnt={}
for i in directory:
if not os.path.exists(dest_dir+i):
os.mkdir(dest_dir+i+'/')
cnt[i]=0

for name in list:
img=Image.open(source_dir+name)
img.crop((5,2,17,20)).save(dest_dir+name[0]+'/'+str(cnt[name[0]]),'png')
cnt[name[0]]+=1
img.crop((17,2,29,20)).save(dest_dir+name[1]+'/'+str(cnt[name[1]]),'png')
cnt[name[1]]+=1
img.crop((29,2,41,20)).save(dest_dir+name[2]+'/'+str(cnt[name[2]]),'png')
cnt[name[2]]+=1
img.crop((41,2,53,20)).save(dest_dir+name[3]+'/'+str(cnt[name[3]]),'png')
cnt[name[3]]+=1

其实也很简单,就是从recognized/中读取源图,分割成4张子图保存在number/中对应的文件夹下,并以数字序数命名。这里用了PIL进行图片分割。但其实这里的难点不是代码本身,而是要观察每张图片分隔的界限。由于普通图片查看器的并没有标尺,为了观察的更细致,我用的是GIMP(Ubuntu中的ps)进行观察。

为了方便纠正训练错误,在每次训练前我会把之前的图片删除重新写入。

package.py

这个文件用来将之前分割过的图片进行转化,也是最关键的部分。

我的思路就是将每个图片用PIL转化成灰度图,在转化为二值图,成为一个12*18的矩阵,再转化为单行数组。数组中每一个点就是一个特征。但是在后来的评估中,为了提高识别准确率,又加入了三个特征,即:空白的联通块数、字符的边界像素个数、字符的填充像素个数。事实证明这三个特征极大的提高了验证的准确性。

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
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import cPickle,os,Image,gzip
import numpy as np
from util import *

print '... packaging'

source_dir='number/'

os.chdir(source_dir)
list=os.listdir('./')
data=[]
ans=[]

dict='012345678abcdefghijklmnpqrstuvwxy'
for num in list:
os.chdir(num)
list_pic=os.listdir('./')
for pic in list_pic:
img=Image.open(pic)
img=blur1(img)
cnt1=seed_fill(img)
cnt2=count_fill(img)
cnt2=(cnt2-50)/10.0
cnt3=count_border(img)
cnt3=(cnt3-50)/10.0
img=blur2(img)
arr=np.array(img)
arr=arr.reshape(12*18).astype('int').astype('bool').astype('float32')
arr=arr.tolist()
arr.append(cnt1)
arr.append(cnt2)
arr.append(cnt3)
arr=np.array(arr).astype('float32')
data.append(arr)
ans.append(float(dict.find(num)))
os.chdir('../')

os.chdir('../')
data=np.array(data)
ans=np.array(ans)

size=len(data)

rand_arr=np.arange(size)
np.random.shuffle(rand_arr)

rand_data=[]
rand_ans=[]

for i in rand_arr:
rand_data.append(data[i])
rand_ans.append(ans[i])

rand_data=np.array(rand_data)
rand_ans=np.array(rand_ans)

len_train=size*0.6
len_valid=size*0.2
len_test=size-len_train-len_valid

output=((rand_data[0:len_train],rand_ans[0:len_train]),(rand_data[len_train:len_train+len_valid],rand_ans[len_train:len_train+len_valid]),(rand_data[size-len_test:size],rand_ans[size-len_test:size]))

outputfile=open('vericode.pkl','wb')
cPickle.dump(output,outputfile)
outputfile.close()

outputfile=open('vericode.pkl','rb')
outputgzipfile=gzip.open('vericode.pkl.gz','wb')
outputgzipfile.writelines(outputfile)
outputfile.close()
outputgzipfile.close()
os.remove('vericode.pkl')

没写注释,就列下注意点吧:

  • 在进行图像格式处理的过程中,很多小细节需要注意,特别是图片到数组的格式转化。
  • 对于新加特征值的计算,我放在了util.py中实现。
  • 对于新加的特征值还要进行特征值的优化,用以提高迭代效率。
  • 对数据需要进行训练集、评估集和测试集分类,比例为6:2:2,并且首先需要进行随机化处理。
  • 输出的验证集我是用cPickle直接输出之后再用gzip进行了压缩。
  • 尤其注意输出的格式,我是用的mnist数据集的标准格式进行处理。

自制正方软件系统验证码的识别程序(1/4)

Posted on 2016-02-05 | In Computer Vision

学了几天机器学习的原理,感觉还是要自己动手写一个像样的程序才行。正好刚研究过deeplearning.net上的那个识别mnist数据库的源码,于是打算利用一下写一个识别验证码的程序。

他提供的代码已经基本实现了利用Logistic Regression 来分类,用的是softmax算法,支持多元分类。然而他只能读入指定的mnist数据集的经过处理后的格式。因此我当前的任务就是将一张实实在在的图片转换成计算机可以识别的格式,并且进行好预处理,选择好特征值,而不用考虑机器学习代码的细节(这个的确有点难)。但是仅仅是处理好这些事情也并不容易。

目标

我的目标是我们学校的选课系统的网站(没错我是苏大的0.0),这里的验证码相对比较规则,而且未来应该有点实用价值吧。

1

就是这里的验证码,取出来看一下就是这样:

当然,为了识别这个,还是得多下载几张下来自己人脑识别几张作为训练集。。。

下载图片

下载这种东西还是用shell脚本方便:

1
2
3
4
for((i=1;i<=10000;i++))
do
wget -O './png/'$i 'http://xk.suda.edu.cn/CheckCode.aspx'
done

一不小心下了10000张,先留着再说。

获得训练集

本来想在网上找点靠谱的识别平台,但是不管要不要钱,效果都不好。所以,,甭想了,先人脑识别500张再说。。。。。。

为了方便以后使用,我就用他的文件名来记录实际的值。

分析

自己观察下载的图片,发现他有这么几个特点:

1、包含字母数字,但是9、o、z三个字符不存在(估计是因为容易产生混淆,为了方便用户就没提供这些)。这一点很重要,因为在进行训练的时候,如果某个集合是空的,那么就会在执行的过程中出现类似‘除0错误’这样的东西。因此分类的时候需要把这三个字符排除。

2、包含四个字符,每一个字符都占据一个相对固定的位置。这一点大大降低了识别的难度,因为这样我们就可以将字符进行分割,从而把问题转化为对单个字符的识别。否则用普通的办法就难以下手了。

3、字符的颜色固定,都是蓝色。这样我们就可以非常方便的进行二值化处理而不用考虑反色的问题了。

4、图片的噪声只有一些象征性的点点。。。这对我们来说几乎形同虚设,只要稍微对图片进行处理就能消除影响。

5、下载了大量的数据之后发现竟然有大量的图片出现重复,于是我猜他的图片并不是动态生成的,而是静态读取的。这么说来就算再不济,我也完全可以把他的图片全部下载,找个验证码平台识别好然后进行hash查找。总共也就167万种,全部下下来就几个G,跑个小半天应该就下的差不多了吧。

根据上面这些特点,我们就能有针对性的设计程序了。

机器学习中的交叉验证思想

Posted on 2016-02-02 | In Machine Learning

简述

在使用训练集对参数进行训练的时候,经常会发现人们通常会将一整个训练集分为三个部分(比如mnist手写训练集)。一般分为:训练集(train_set),评估集(valid_set),测试集(test_set)这三个部分。这其实是为了保证训练效果而特意设置的。其中测试集很好理解,其实就是完全不参与训练的数据,仅仅用来观测测试效果的数据。而训练集和评估集则牵涉到下面的知识了。

因为在实际的训练中,训练的结果对于训练集的拟合程度通常还是挺好的(初试条件敏感),但是对于训练集之外的数据的拟合程度通常就不那么令人满意了。因此我们通常并不会把所有的数据集都拿来训练,而是分出一部分来(这一部分不参加训练)对训练集生成的参数进行测试,相对客观的判断这些参数对训练集之外的数据的符合程度。这种思想就称为交叉验证(Cross Validation)。

通常我们使用的交叉验证方法有下面几种:

简单交叉验证(simple cross validation)

简单交叉验证当然很简单了,就是把整个训练集随机分为两部分(通常是70%的训练集,30%的评估集)。

1、首先我们用训练集建立模型,这时候我们需要尝试多种的参数来得到一些不同的模型;

2、对于每一个模型,调用评估集来进行测试,计算出训练误差(通常是以类似损失函数的形式);

3、取训练误差最小的那个模型作为最后的结果;

这个方法当然会存在一些训练数据的浪费,毕竟有大量的数据没有参与真正的训练而是仅仅作为评估。所以这个方法只能在数据非常易得的情况下使用,如果数据比较珍贵,显然这种方法就不适用了。

有时候这个方法好像也被称为HoldOut验证(Hold-Out Method)。其实这也不算是交叉验证了,因为他的训练集并没有交叉。

通常情况下我们是直接选取前70%为训练集,但是如果训练数据是按照一定规律排放的,那么选取数据的时候就要先打乱顺序,或者按照一定的随机方法选取数据。否则训练集就不一定具有一般性了。

K-折交叉验证(S-fold Cross Validation)

这个据说是最常用的验证方法了,步骤如下:

1、将数据集均分为K份

2、从K份中取一份作为评估集,另外K-1份作为训练集,生成K个模型以及这K个模型对于评估集的训练误差;

3、取训练误差最小的那个模型作为最后的结果;

经大量实验验证,据说我们取K=10的时候效果最好。

这个方法一方面保证了数据充分被使用训练了,避免了数据的浪费;另一方面也互相进行了验证,达到了交叉验证的效果,不过计算代价还是有点高。

留p交叉验证(Leave-p-out Cross Validation)

从名字大概就可以看出来了,所谓留p,就是每一次训练都会留下p个数据作为评估集,剩下的n-p个数据作为训练集,分别进行建模测试,取出效果最好的模型。这里我们会发现这p个数据可以有$C_n^p$组可能,所以计算的复杂度还是挺高的。

当然,有时候我们会取p=1,这导致了每一个数据都会独立作为测试集。这种方法又被叫做**留一交叉验证(Leave-One-Out Cross Validation)**,当数据极为匮乏的时候才会使用。

事实上,交叉验证的方法不仅能够提高数据的利用率,更重要的是他也能够在一定程度上解决过拟合(Overfitting)问题,因为过拟合只能很好的拟合训练集中的数据而并不能拟合评估集中的数据。

1…464748…58

574 posts
69 categories
286 tags
© 2024 Companyd
Powered by Hexo
|
Theme — NexT.Muse v5.1.4