% 數字識別
% 引自《神經網絡模型及其MATLAB仿真程序設計》一書P94～98

% 生成輸入向量和目標向量
clear all;
'LOADING......'

for kk = 0:99
    p1=ones(16,16);
    m=strcat('nums\',int2str(kk),'.bmp');
    x=imread(m,'bmp');
    bw=im2bw(x,0.5);
    [i,j]=find(bw==0);
    imin=min(i);
    imax=max(i);
    jmin=min(j);
    jmax=max(j);
    bw1=bw(imin:imax,jmin:jmax);
    rate=16/max(size(bw1));
    bw1=imresize(bw1,rate);
    [i,j]=size(bw1);
    i1=round((16-i)/2);
    j1=round((16-j)/2);
    p1(i1+1:i1+i,j1+1:j1+j)=bw1;
    p1=-1.*p1+ones(16,16);
    for m=0:15
        p(m*16+1:(m+1)*16,kk+1)=p1(1:16,m+1);
    end
    
    switch kk
        case{0,10,20,30,40,50,60,70,80,90}
            t(kk+1)=0;
        case{1,11,21,31,41,51,61,71,81,91}
            t(kk+1)=1;
        case{2,12,22,32,42,52,62,72,82,92}
            t(kk+1)=2;
        case{3,13,23,33,43,53,63,73,83,93}
            t(kk+1)=3;
        case{4,14,24,34,44,54,64,74,84,94}
            t(kk+1)=4;
        case{5,15,25,35,45,55,65,75,85,95}
            t(kk+1)=5;
        case{6,16,26,36,46,56,66,76,86,96}
            t(kk+1)=6;
        case{7,17,27,37,47,57,67,77,87,97}
            t(kk+1)=7;
        case{8,18,28,38,48,58,68,78,88,98}
            t(kk+1)=8;
        case{9,19,29,39,49,59,69,79,89,99}
            t(kk+1)=9;
    end
end

'LOAD OK.'

save E52PT p t;
% 創建和訓練BP網絡

clear all;
load E52PT p t;

pr(1:256,1)=0;
pr(1:256,2)=1;
net=newff(pr,[25 1],{'logsig' 'purelin'}, 'traingdx', 'learngdm');
net.trainParam.epochs=2500;
net.trainParam.goal=0.001;
net.trainParam.show=10;
net.trainParam.lr=0.05;
net=train(net,p,t)

'TRAIN OK.'

save E52net net;

% 識別
for times=0:999
    clear all;
    p(1:256,1)=1;
    p1=ones(16,16);
    load E52net net;
    test=input('FileName:', 's');
    x=imread(test,'bmp');
    bw=im2bw(x,0.5);
    [i,j]=find(bw==0);
    imin=min(i);
    imax=max(i);
    jmin=min(j);
    jmax=max(j);
    bw1=bw(imin:imax,jmin:jmax);
    rate=16/max(size(bw1));
    bw1=imresize(bw1,rate);
    [i,j]=size(bw1);
    i1=round((16-i)/2);
    j1=round((16-j)/2);
    p1(i1+1:i1+i,j1+1:j1+j)=bw1;
    p1=-1.*p1+ones(16,16);
    for m=0:15
        p(m*16+1:(m+1)*16,1)=p1(1:16,m+1);
    end
    [a,Pf,Af]=sim(net,p);
    imshow(p1);
    a=round(a)
end