目前，pca算法已经广泛应用于各方面，就拿图像处理，经常做的一件事就是当提取的图像特征维度比较高时，为了简化计算量以及储存空间，需要对这些高维数据进行一定程度上的降维，并尽量保证数据的不失真。
先举个例子，方便理解：
1）对于一个训练集，20个sample(i=1,2,3,...,20)，特征Xi是100维[Xi1,Xi2,Xi3,...Xij,...,Xi100](j=1,2,..,100)，那么它可以建立一个20*100的样本矩阵M。
2）紧接着我们开始求这个样本的协方差矩阵，得到一个20*20的协方差矩阵，计算过程如下：
•先求解出Xi的平均Xav=(∑xi)/20;
•对每一个Xi,计算Xi-Xav,即Mi(第 i 行)变为 Mi-Xav,记为Mn;
•则容易得到协方差矩阵Z为Mn*Mn'( ' 表示转置 ) 。
3）然后求出这个协方差矩阵Z20x20的特征值和特征向量，一般情况下应该有20个特征值和特征向量，现在根据特征值的大小，取出较大的特征值以及其所对应的特征向量，（假设提取的特征值为较大的5个特征值），那么这5个特征向量就会构成一个20*5的矩阵V，这个矩阵就是我们要求的特征矩阵。
4）用Mn'去乘以V，得到一个base矩阵（*），大小为100x5。
5）任取一个样本1x100,乘上这个100*5的特征矩阵，就得到了一个1*5的新的样本，显然每个sample的维数下降了，然后再用这个1x5向量去比较相似性。
注：
›上述3）过程中特征值的选取在不确定具体要降到多少维的情况下，一般还可以根据n个特征值之和大于总和的90%进行选取。
›上面的（*）处base矩阵的求解不唯一，也可以自行修正。
大致说完了PCA降维的过程，现在会有人要问为什么只提取特征值较大的几个特征值就可以近似代替原样本矩阵呢。
好了，话不多说。下面就讲讲矩阵的特征值和特征向量的数学意义：
为简单起见，以二维矩阵A=[1 0;0 -1]（矩阵的秩为2）为例，以平面上一点（x,y）经过A变换后变为（x',y'）若这两点在一条直线在，那么可以理解为矩阵A的作用恰好使得向量[x y]' 只是在原有方向上变换了长度而已，即Ax=λx (x为一列向量).对于A矩阵，容易得到A的两个特征值及相应的特征向量 λ1=1 ,e1=[1 0]' , λ2=-1 ,e2=[0 -1]' ,二维平面上任意一点（x,y）=b1*(1,0)+b2*(0,-1)（b1,b2均为实常数）; 那么A[x y]'=A*(b1*e1+b2*e2)=b1*λ1+b2*λ2 =∑biλi ;
把这个公式推广到高维空间，在计算（x',y'）对于λ值比较小的特征维可以忽略.
B=[1 0;0 0.01] ,其中B的两个特征值及相应的特征向量 λ1=1 ,e1=[1 0]' , λ2=0.01 ,e2=[0 1]'
那么x=[2 3]' 经过B变换为 Bx=[2 0.03]';
如果我们认为λ2远小于λ1，忽略掉λ2时，Bn=[1 0;0 0],Bnx=[2 0]'≈[2 0.03].
通俗点讲，pca算法就是去寻找那些在该维度上方差比较大的维，同时忽略比较平均的维度。假如上面所说的X特征向量的第一个元素都为1，那么这一列数据是可以忽略的，因为他并不能起到区分的作用，相反我们是要寻找那些在某一维度分布比较广的维，并提取出来。
打个比方，平面区域一个斜75度的椭圆，且长轴远大于短轴，那么椭圆上的点在短轴上的分布明显弱于长轴，当短轴远小于长轴时，近似为一条直线，便失去了短轴这个维度。
下面贴上一个matlab版的pca算法代码：
%一个修改后的PCA进行人脸识别的Matlab代码
clear;
% calc xmean,sigma and its eigen decomposition
allsamples=[];%所有训练图像
for i=1:40
for j=1:5
a=imread(strcat('C:\Documents and Settings\Foreigners\桌面\ORL\s',num2str(i),'\',num2str(j),'.bmp'));
% imshow(a);
b=a(1:112*92); % b是行矢量 1×N，其中N＝10304，提取顺序是先列后行，即从上到下，从左到右
b=double(b);
allsamples=[allsamples; b]; % allsamples 是一个M * N 矩阵，allsamples 中每一行数据代表一张图片，其中M＝200
end
end
samplemean=mean(allsamples); % 平均图片，1 × N
for i=1:200 xmean(i,:)=allsamples(i,:)-samplemean; % xmean是一个M × N矩阵，xmean每一行保存的数据是“每个图片数据-平均图片”
end;
sigma=xmean*xmean'; % M * M 阶矩阵
[v d]=eig(sigma);
d1=diag(d);
[d2 index]=sort(d1); %以升序排序
cols=size(v,2);% 特征向量矩阵的列数
for i=1:cols
vsort(:,i) = v(:, index(cols-i+1) ); % vsort 是一个M*col(注:col一般等于M)阶矩阵，保存的是按降序排列的特征向量,每一列构成一个特征向量
dsort(i) = d1( index(cols-i+1) ); % dsort 保存的是按降序排列的特征值，是一维行向量
end %完成降序排列
%以下选择90%的能量
dsum = sum(dsort);
dsum_extract = 0;
p = 0;
while( dsum_extract/dsum < 0.90)
p = p + 1;
dsum_extract = sum(dsort(1:p));
end
i=1;
% (训练阶段)计算特征脸形成的坐标系
while (i<=p && dsort(i)>0)
base(:,i) = dsort(i)^(-1/2) * xmean' * vsort(:,i); % base是N×p阶矩阵，除以dsort(i)^(1/2)是对人脸图像的标准化，详见《基于PCA的人脸识别算法研究》p31
i = i + 1;
end
% add by wolfsky 就是下面两行代码，将训练样本对坐标系上进行投影,得到一个 M*p 阶矩阵allcoor
allcoor = allsamples * base;
accu = 0;
% 测试过程
for i=1:40
for j=6:10 %读入40 x 5 副测试图像
a=imread(strcat('C:\Documents and Settings\Foreigners\桌面\ORL\s',num2str(i),'\',num2str(j),'.bmp'));
b=a(1:10304);
b=double(b);
tcoor= b * base; %计算坐标，是1×p阶矩阵
for k=1:200
mdist(k)=norm(tcoor-allcoor(k,:));
end;
%三阶近邻
[dist,index2]=sort(mdist);
class1=floor( index2(1)/5 )+1;
class2=floor(index2(2)/5)+1;
class3=floor(index2(3)/5)+1;
%class=class1;%%blue_lg
if class1~=class2 && class2~=class3
class="class1";
elseif class1==class2
class="class1";
%elseif class2==class3
% class="class2";
end;
if class==i
accu=accu+1;
end;
end;
end;
accuracy=accu/200 %输出识别率
%zuobiao=[1:100];
%plot(zuobiao,accuracy);
转载请注明： http://www.cnblogs.com/blue-lg/archive/2012/05/14/2499581.html
下面介绍一个自己搞的小时钟，方便拖拽。
那么如何插入到自己的界面中呢？
只需要加上如下代码就好了：
<script type="text/javascript" src="clock.js"> </script>
clcok.js的完整代码如下：
var extra=document.createElement('div');
extra.style.position='absolute';
var extra_canvas=document.createElement('canvas');
extra_canvas.id="extra_canvas";
extra.appendChild(extra_canvas);
document.body.appendChild(extra);
var flag;
var currentRectX;
var currentRectY;
function init(){
flag=true;
currentRectX=0;
currentRectY=0;
}
function clock(size,x,y){
/* if(!flag){
document.body.removeChild(extra_canvas);
console.log('remove');
}
*/
if (!size){size=200;}
if (!x)x=0;
if (!y)y=0;
extra_canvas.width=size;
extra_canvas.height=size;
var context=extra_canvas.getContext('2d');
extra.style.left=currentRectX+'px';
extra.style.top=currentRectY+'px';
//console.log(currentRectX,currentRectY);
//context.fillStyle='#FF0000';
//context.fillRect(0,0,100,100);
//绘制表盘
var centerX=x+size/2;
var centerY=y+size/2;
//console.log(centerX,centerY);
var radius=(size-40)/2;
context.clearRect(x,y,x+size,y+size);
context.beginPath();
context.arc(centerX, centerY, radius, 0, 2 * Math.PI, false);
context.lineWidth = 5;
context.strokeStyle = "grey";
context.stroke();
var grd = context.createLinearGradient(centerX-radius, centerY-radius,centerX+radius, centerY+radius);
grd.addColorStop(0, "#FFFFFF"); // light blue
grd.addColorStop(1, "#DDDDFF"); // dark blue
context.fillStyle = grd;
context.fill();
context.closePath();
//绘制数字时刻
context.beginPath();
//context.font="9pt";
// context.fontsize=9px;
//context.fontFamily="Square721 BT";
context.fillStyle = "#0000f0"; // text color
context.fillText("12",centerX-7,centerY-radius+18);
context.fillText("3",centerX+radius-18,centerY+4);
context.fillText("6",centerX-3,centerY+radius-12);
context.fillText("9",centerX-radius+12,centerY+4);
context.closePath();
/* //显示日期
context.beginPath();
context.TextOut(100,100,"星期");
context.closePath();
*/
//绘制刻度
for (var i=0;i<12;i++){
context.beginPath();
if(i%3){
context.lineWidth = 3;
context.strokeStyle = "grey";
len=5;
}else{
context.lineWidth = 6;
context.strokeStyle = "#ff0";
len=10;
}
arc=i/6*Math.PI;
kx=centerX+radius*Math.cos(arc);
ky=centerY-radius*Math.sin(arc);
context.moveTo(kx,ky);
kx=centerX+(radius-len)*Math.cos(arc);
ky=centerY-(radius-len)*Math.sin(arc);
context.lineTo(kx,ky);
context.stroke();
context.closePath();
}
//绘制表中心
context.beginPath();
context.arc(centerX, centerY, 4, 0, 2 * Math.PI, false);
context.lineWidth = 2;
context.strokeStyle = "grey";
context.stroke();
context.closePath();
//绘制时针分针
var now=new Date();
var hour=now.getHours()%12;
var minute=now.getMinutes();
var second=now.getSeconds();
hour=hour+minute/60;//update the time!!
minute=minute+second/60;
var arc_hour=hour/6*Math.PI;
context.beginPath();
kx=centerX+(radius-40)*Math.sin(arc_hour);
ky=centerY-(radius-40)*Math.cos(arc_hour);
context.moveTo(kx,ky);
kx=centerX+10*Math.sin(arc_hour+Math.PI);
ky=centerY-10*Math.cos(arc_hour+Math.PI);
context.lineTo(kx,ky);
context.lineWidth = 6;
context.strokeStyle = "black";
context.stroke();
context.closePath();
context.beginPath();
var arc_minute=minute/30*Math.PI;
context.beginPath();
kx=centerX+(radius-20)*Math.sin(arc_minute);
ky=centerY-(radius-20)*Math.cos(arc_minute);
context.moveTo(kx,ky);
kx=centerX+20*Math.sin(arc_minute+Math.PI);
ky=centerY-20*Math.cos(arc_minute+Math.PI);
context.lineTo(kx,ky);
context.lineWidth = 4;
context.strokeStyle = "red";
context.stroke();
context.closePath();
// flag=false;
context.beginPath();
var arc_second=second/30*Math.PI;
context.beginPath();
kx=centerX+(radius-20)*Math.sin(arc_second);
ky=centerY-(radius-20)*Math.cos(arc_second);
context.moveTo(kx,ky);
kx=centerX+20*Math.sin(arc_second+Math.PI);
ky=centerY-20*Math.cos(arc_second+Math.PI);
context.lineTo(kx,ky);
context.lineWidth = 2;
context.strokeStyle = "gray";
context.stroke();
context.closePath();
}
/*
extra.onmousedown=null;
extra.onmouseup=null;
extra.onmousemove=null;
*/
extra_canvas.onmousedown=canvasMouseDownHandler;
extra_canvas.onmouseup=canvasMouseUpHandler;
function canvasMouseDownHandler(event){
var canvasMouseX=event.screenX;
var canvasMouseY=event.screenY;
extra_canvas.onmousemove=canvasMouseMoveHandler;
// console.log('down');
startDragMouseX=canvasMouseX;
startDragMouseY=canvasMouseY;
startDragRectX=currentRectX;
startDragRectY=currentRectY;
}
function canvasMouseMoveHandler(event){
event.preventDefault();
var canvasMouseX=event.screenX;
var canvasMouseY=event.screenY;
// console.log('move');
currentRectX=startDragRectX+canvasMouseX-startDragMouseX;
currentRectY=startDragRectY+canvasMouseY-startDragMouseY;
//console.log(currentRectX,currentRectY);
clock();
}
function canvasMouseUpHandler(event){
extra_canvas.onmousemove=null;
//console.log('up');
}
/* function cc(){
clock(null,0,0);
}
*/
//window.setInterval(cc, 200);
init();
window.setInterval(clock, 200);
代码部分的注释大家可以忽略哈~
效果嘛，看http://blue-lg.diandian.com/页面~~~
谢谢大家关注，转载请注明：http://www.cnblogs.com/blue-lg/archive/2012/03/22/2411490.html
时至今日，jpg图像文件成为了一般网络传递的首选图片格式。相比于bmp文件，jpg像素是有损失的。
下面我们开始了解jpeg文件的压缩编码过程。
首先我们得了解DCT离散余弦变换和霍夫曼编码以及其他一些常用的编码方式。
JPEG的压缩编码过程如下：
1.颜色模型转换
这一步是将常见的RGB颜色模型转变为YCrCb模型，
转换公式如下：
Y=0.299R+0.587G+0.114B;
Cb=-0.1687R-0.3313G+0.5B+128;
Cr=0.5R-0.4187G-0.0813B+128.
[反推：R=Y+1.402(Cr-128);G=Y-0.34414(Cb-128)-0.71414(Cr-128);B=Y+1.772(Cb-128).]
至于原因,看了下面的例子，想必能有点启发。
由于人眼对亮度Y远比色度C更敏感，所以如果2X2个点表示成RGB格式，需要4*3*1=12个字节，如果用4:2:0YCrCb模型，这几个点的Cr和Cb近似认为相等，即4个点公用一对色度，这样占用的字节就变为4+2=6字节，这大大减小了存储空间。
关于4:2:0模型，相关的还有4:2:1等，详情查看维基百科-色度抽样。
2.分块并进行DCT变换
将M*N的图像矩阵分为8*8块存储，基于图像像素的渐变性，可以利用差分编码减小字节数。
对每一块8*8数组矩阵进行DCT离散余弦变换，
其中u、v分别表示像素点的坐标位置，当u=v=0时,C(u)=C(v)=1/1.414,反之，C(u)=C(v)=1。（N=8）经过上面的公式变换，得到8*8=64个系数，一般而言这些系数中以F(0,0)的值最大，显而易见，以亮度C为例，F(0,0)代表的是原本块的平均亮度，即直流分量DC。剩下的63个系数多半是接近于0的值，他们称为AC分量。
但是，由于单纯的进行一次变换并不能使得系数的个数变少，始终为64，但是不同的是，有大多数值集中在横轴0上下，所以下一步操作就是量化，将这些微小的高频系数“忽略”掉。
3.量化
所谓量化，指的就是将系数按比例缩小，比如下面一组数据A=[1，14，7，19，16，23，12]
均以4作为量化步长，得到的新的一组数据为B=[0，3，2，4，6，3]，解码的时候反量化得到A'=[0，12，8，16，24，12].
量化器的输出值取最接近被量化步长整除的的整数作为新的数据。
JPEG编码使用的量化指表如下所示：
亮度：
注意C和Y采用的是不同的量化表，这里没有介绍色度的量化表，读者可自行查询。
接下来，就开始着手准备数据的编码，这里采用Z编码。
4.编码
如上图所示，将8*8数组按Z顺序排布读取系数，这样读取可以增加值连续为0的长度，减小储存空间。
对于直流系数和交流系数分别采用不同的编码方式，
直流系数：
利用相邻块DC值差别较小，使用差分编码，对于每个块的DC值取差值，D'=DC(i)-DC(i-1);
交流系数：
使用游程编码方法，对上图的系数进行游程编码，得到结果如下：
理想化系数：-3,0,5,0,0,0,9,0,4,0,0,...
编码结果：（0，-3），（1，5），（3，9），（1，4），（0，0）。
其中（0,0）代表结束标志，（a，b）表示当前值与前一个非零值之间的0的个数，b则为当前非零值。
将上面得到的亮度、色度的DC和AC4种码统一利用霍夫曼编码进行最终的编码。
读者可以在网络上查询上述的4种不同的霍夫曼表，用来对系数进行编码。（图像编码基础和小波压缩技术 第二章附录）
解码过程需要注意的是jpeg文件读取的几个重要标识。
1. 读取文件的大体结构
JFIF格式的JPEG文件(*.jpg)的一般顺序为：
SOI(0xFFD8),APP0(0xFFE0),[APPn(0xFFEn)]可选,
DQT(0xFFDB),SOF0(0xFFC0),DHT(0xFFC4),SOS(0xFFDA),
压缩数据,EOI(0xFFD9)。
2. 读取哈夫曼表数据；
3. 建立哈夫曼树。
在准备好所有的图片信息后，就可以对图片数据进行解码了。
过程就是把上述几个步骤反过来弄一遍。
图像的细化是模式识别中很重要的一个技术，指的是将原本"臃肿"的像素简化为单像素相连接的二值图像（即类似骨架的概念），细化的好坏直接影响到后面识别匹配的效率。
摘自某文章的话，细化就是经过一层层的剥离，从原来的图中去掉一些点，但仍要保持原来的形状，直到得到图像的骨架。骨架，可以理解为图象的中轴，例如一个长方形的骨架是它的长方向上的中轴线；正方形的骨架是它的中心点；圆的骨架是它的圆心，直线的骨架是它自身，孤立点的骨架也是自身。
下面先介绍经典的Zhang并行快速细化算法：
设p1点的八邻域为：
【 p9 p2 p3
p8 p1 p4
p7 p6 p5 】
（其中p1为白点，如果以下四个条件同时满足，则删除p1，即令p1=0）
其中迭代分为两个子过程：
过程1 细化删除条件为： （1）、2 < =N(p1) <= 6， N(x)为x的8邻域中黑点的数目
（2）、A(p1)=1, A(x)指的是将p2-p8之间按序前后分别为0、1的对数（背景色：0）
（3）、p2*p4*p6=0
（4）、p4*p6*p8=0
如果同时满足以上四个条件则该点可以删除（赋值为0）。
过程2 细化删除条件为： （1）、2 < =N(p1) <= 6， N(x)为x的8邻域中黑点的数目
（2）、A(p1)=1, A(x)指的是将p2-p8之间按序前后分别为0、1的对数（背景色：0）
（3）、p2*p4*p8=0
（4）、p2*p6*p8=0
如果同时满足以上四个条件则该点可以删除。
代码如下：
A.m
1 function n=A(temp,i,j)
2 %0->1的数目
3 shuzu=[temp(i,j),temp(i-1,j),temp(i-1,j+1),temp(i,j+1),temp(i+1,j+1),temp(i+1,j),temp(i+1,j-1),temp(i,j-1),temp(i-1,j-1)];
4 n=0;
5 for i=2:8
6 if shuzu(i)==0&&shuzu(i+1)==1
7 n=n+1;
8 end
9 end
主函数代码：
1 test=input('Please input a digits image:','s'); %输入图像
2 x=imread(test);
3 if ~isbw(x)
4 '请确保输入图像为二值化图像！';
5 else
6 [height,width]=size(x);
7 mark=1;
8 % temp=zeros(height+2,width+2);
9 % temp(2:height+1,2:width+1)=x(:,:);
10 temp=x;
11 imshow(temp);
12 while mark==1
13 mark=0;
14
15 for i=2:height-1
16 for j=2:width-1
17 condition=0;
18 %判断P(r,c)是否为可细化像素
19 if temp(i,j)==1
20 n=0;
21 for ii=-1:1
22 for jj=-1:1
23 n=n+temp(i+ii,j+jj);
24 end
25 end
26 if (n>=3 && n<=7)
27 condition=condition+1;
28 end
29 if A(temp,i,j)==1
30 condition=condition+1;
31 end
32 if temp(i-1,j)*temp(i,j+1)*temp(i+1,j)==0
33 condition=condition+1;
34 end
35 if temp(i,j+1)*temp(i+1,j)*temp(i,j-1)==0
36 condition=condition+1;
37 end
38 if condition==4
39 mark=1;
40 temp(i,j)=0;
41 end
42 end
43 end
44 end
45 figure;imshow(temp);
46
47
48 for i=2:height-1
49 for j=2:width-1
50 condition=0;
51 %判断P(r,c)是否为可细化像素
52 if temp(i,j)==1
53 n=0;
54 for ii=-1:1
55 for jj=-1:1
56 n=n+temp(i+ii,j+jj);
57 end
58 end
59 if (n>=3 && n<=7)
60 condition=condition+1;
61 end
62 if A(temp,i,j)==1
63 condition=condition+1;
64 end
65 if temp(i-1,j)*temp(i,j+1)*temp(i,j-1)==0
66 condition=condition+1;
67 end
68 if temp(i,j-1)*temp(i+1,j)*temp(i,j-1)==0
69 condition=condition+1;
70 end
71 if condition==4
72 mark=1;
73 temp(i,j)=0;
74 end
75 end
76 end
77 end
78 figure;imshow(temp);
79 end
80 end
结果：
离线手写数字识别系统设计主要分为3个部分：图像预处理、数字的分割、数字的识别。
1.图像预处理
图像预处理主要是指对输入的图像进行灰度化、滤波、二值化。
鉴于数字的识别与色彩无关，并且考虑到噪声的影响（这里采用中值滤波），将图像进行预处理。
1 clear all;
2 clc;
3 test=input('Please input a digits image:','s'); %输入图像
4 x=imread(test);
5 if ~isgray(x)
6 x=rgb2gray(x); %必须转换为灰度图像
7 end
8 xbw=im2bw(x,0.9); %再转换为二值图像
9 xbw=medfilt2(xbw); %中值滤波
10 bw=xbw; %滤波后二值图像
11
12 figure;imshow(bw);
2.图像切割（字符分离）
原理：
当图像转换为二值图像后，利用边界点确定单个字符（假设每个数字之间均存在空白、间隔），如下：
当需要提取‘4’这个字符的时候，需要在横向、纵向上确定字符的上下边界点。
这里需要用到函数find()和all()【find()结果为满足条件的所有值，all()：所有的点均满足某条件的值】
先确定左边界，再找到全1列（右边红线），以此类推确定上下边界。
1 for n=1:100
2 [c,l]=size(bw);
3 [i,j] = find(bw==0);
4 jmin1=min(j);
5 xbwt=bw(:,jmin1:l); %xbwt为切除左边边缘的图像
6 [c,l]=size(xbwt);
7 [i,j]=find(all(xbwt(1:c,:)==1));
8 jmin2=min(j); %找右边界
9 if(isempty(jmin1)||isempty(jmin2)) break;
10 end
11 xbwn{n}=xbwt(:,1:jmin2-1); %xbwn{n+1}为第（n+1）个数字图像
12 bw=xbwt(:,jmin2:l); %bw迭代循环使用
13 [i,j] = find(xbwn{n}==0); %找上下边界
14 imin=min(i);
15 imax=max(i); %imin:数字最上点,imax:数字最下点
16 xbwn1{n}=xbwn{n}(imin:imax,:);
17
18 end
如上图所示数字太过饱满，需要骨架化一下才能方便识别，这里需要用到matlab的bwmorph()函数，详见http://www.mathworks.cn/help/toolbox/images/ref/bwmorph.html
1 m=n-1;
2 for n=2:m
3 [i,j]=size(xbwn1{n});
4 p=-1.*xbwn1{n}+ones(i,j); %先将矩阵取反
5 xbwn2{n}=bwmorph(p,'skel',Inf);
6 %figure,imshow(xbwn2{n});
7 end
3.字符的匹配识别
1）统一标准
将提取出来的字符图片大小全部统一为16X16
2）提取特征向量
特征向量x为由以下数据组成的数组，
x1：数字的端点数（8邻域内黑点数目为3（包括本身）），
x2：数字一般连接点（字轨迹长度），
x3：数字分叉点（8邻域内黑点数目大于3（包括本身）），
x4~x19：16条水平线分别与数字的交点个数，
x20~x35：16条垂直线分别与数字的交点个数 。
3）匹配确定数字
假设标准数字(0~9)的特征向量已知为y，则匹配的准则为计算x和y的相似度r.
r=Σ(xi*yi)/sqrt(Σ(xi^2)*Σ(yi^2));【sqrt为平方根运算】
r数字越接近于1，表明y和x越相似，即未知字符最有可能为y对应的标准数字。
下面是matlab代码:
solve.m
1 function x=solve(ju)
2 x=zeros(1,35);
3 [m,n]=size(ju);
4 ju1=zeros(m+2,n+2);
5 ju1(2:m+1,2:n+1)=ju(:,:);
6 for i=2:m+1
7 for j=2:n+1
8 if ju1(i,j)==1
9 temp=ju1(i-1:i+1,j-1:j+1);%邻域矩阵
10 [ii,jj]=find(temp==1);
11 if max(size(ii))>=4
12 x(3)=x(3)+1;
13 else if max(size(ii))==3
14 x(2)=x(2)+1;
15 else
16 x(1)=x(1)+1;
17 end
18 end
19 x(i+2)=x(i+3)+1;
20 x(j+18)=x(j+19)+1;
21 end
22 end
23 end
值得注意的是，上面第5行用ju1代替ju，是为了方便判断边缘上的点（ju此时寻找8邻域会出错），所以在四周加上了无关‘边’。
pipei.m
1 function result=pipei(a,j)
2 sta=[0 27 5 0 1 3 2 2 2 2 2 2 2 2 2 2 2 2 3 1 0 0 0 0 0 11 5 2 2 3 9 0 0 0 0
3 2 14 0 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 15 1 0 0 0 0 0 0
4 2 18 3 0 0 0 3 2 2 1 1 1 1 1 1 1 1 1 5 2 0 0 0 0 0 6 4 3 3 4 3 0 0 0 0
5 3 18 4 0 0 3 2 1 1 1 1 4 1 1 1 2 2 2 3 0 0 0 0 0 0 1 3 3 3 3 9 3 0 0 0
6 3 20 6 0 1 1 1 1 2 2 2 2 2 2 2 2 6 1 1 1 0 0 0 0 0 4 4 2 4 13 1 1 0 0 0
7 3 23 3 0 5 1 1 1 1 4 2 2 1 1 1 2 2 2 3 0 0 0 0 0 0 6 4 3 3 3 9 1 0 0 0
8 1 24 4 0 0 3 2 1 1 1 5 2 2 2 2 2 2 2 2 0 0 0 0 0 0 11 3 3 3 4 5 0 0 0 0
9 2 17 3 0 6 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 1 1 5 6 4 4 1 0 0 0
10 0 23 13 0 1 4 2 2 2 2 2 4 2 2 2 2 2 2 4 1 0 0 0 0 0 8 8 3 3 6 8 0 0 0 0
11 1 28 4 0 0 4 2 2 2 2 2 2 2 3 3 1 2 2 4 0 0 0 0 0 0 9 3 3 3 3 12 0 0 0 0];
12 sx2=0;sy2=0;sxy=0;
13 for i=1:35
14 sx2=sx2+a(i)*a(i);
15 sy2=sy2+sta(j,i)*sta(j,i);
16 sxy=sxy+a(i)*sta(j,i);
17 end
18 result=sxy/sqrt(sx2*sy2);
19
上面的sta数据为测试所得，由于标准有异，读者可以自己收集数据，调试。
完整的主程序如下：
View Code
1 clear all;
2 clc;
3 test=input('Please input a digits image:','s'); %输入图像
4 %z=input('z:','s');
5 x=imread(test);
6 if ~isgray(x)
7 x=rgb2gray(x); %必须转换为灰度图像
8 end
9 xbw=im2bw(x,0.9); %再转换为二值图像
10 xbw=medfilt2(xbw); %中值滤波
11 bw=xbw; %滤波后二值图像
12
13 figure;imshow(bw);
14
15 for n=1:100
16 [c,l]=size(bw);
17 [i,j] = find(bw==0);
18 jmin1=min(j);
19 xbwt=bw(:,jmin1:l); %xbwt为切除左边边缘的图像
20 [c,l]=size(xbwt);
21 [i,j]=find(all(xbwt(1:c,:)==1));
22 jmin2=min(j); %找右边界
23 if(isempty(jmin1)||isempty(jmin2)) break;
24 end
25 xbwn{n}=xbwt(:,1:jmin2-1); %xbwn{n+1}为第（n+1）个数字图像
26 bw=xbwt(:,jmin2:l); %bw迭代循环使用
27 [i,j] = find(xbwn{n}==0); %找上下边界
28 imin=min(i);
29 imax=max(i); %imin:数字最上点,imax:数字最下点
30 xbwn1{n}=xbwn{n}(imin:imax,:);
31
32 end
33 m=n-1;
34 for n=2:m
35 [i,j]=size(xbwn1{n});
36 p=-1.*xbwn1{n}+ones(i,j); %先将矩阵取反
37 x=16;
38 p1=ones(x,x);
39 %xbwn2{n}=bwmorph(p,'skel',Inf);
40 xbwn2{n}=xbwn1{n};
41 %转换到x*x矩阵
42 rate=x/max(size(xbwn2{n}));
43 xbwn3{n}=imresize(xbwn2{n},rate);
44 [i,j]=size(xbwn3{n});
45 i1=round((x-i)/2);
46 j1=round((x-j)/2);
47 p1(i1+1:i1+i,j1+1:j1+j)=xbwn3{n};
48 p1=ones(x,x)-p1;
49 p1=bwmorph(p1,'skel',Inf);
50 p2{n}=p1;
51
52 temp1=solve(p1);
53 temp2=0;
54 temp_index=0;
55 for tt=1:10
56 temp3=pipei(temp1,tt);
57 if temp3>temp2
58 temp2=temp3;
59 temp_index=tt-1;
60 end
61 end
62 temp_index %显示数字
63
64 figure,imshow(p1);
65 end
参考论文：
基于图像识别技术的手写数字识别方法 吴忠1，2，朱国龙1，2，黄葛峰1，2，吴建国1，2
需要完善的地方：
1.评判机制
可以适当调整x1~x3等的权系数，使得相似度测量更精确；
2.标准手写字体的数据收集
本实验只是找了一个印刷体的数据得到了sta数组；
3.字符的分割
中值滤波不一定能很好的滤去一些端状杂质，需要根据高度和宽度大小限制再进一步滤波。
最后，有不当之处，希望园友指正。
转载请注明：http://www.cnblogs.com/blue-lg/archive/2012/03/05/2381078.html