clear ;
close all;
%Step1 獲取圖像   裝入待處理彩色圖像并顯示原始圖像
Scolor = imread('3.jpg');%imread函數讀取圖像文件
%將彩色圖像轉換為黑白并顯示
Sgray = rgb2gray(Scolor);%rgb2gray轉換成灰度圖
figure,imshow(Scolor),title('原始彩色圖像');%figure命令同時顯示兩幅圖像
figure,imshow(Sgray),title('原始黑白圖像');
%Step2 圖像預處理   對Sgray 原始黑白圖像進行開操作得到圖像背景
s=strel('disk',13);%strei函數
Bgray=imopen(Sgray,s);%打開sgray s圖像
figure,imshow(Bgray);title('背景圖像');%輸出背景圖像
%用原始圖像與背景圖像作減法，增強圖像
Egray=imsubtract(Sgray,Bgray);%兩幅圖相減
figure,imshow(Egray);title('增強黑白圖像');%輸出黑白圖像
%Step3 取得最佳閾值，將圖像二值化
fmax1=double(max(max(Egray)));%egray的最大值并輸出雙精度型
fmin1=double(min(min(Egray)));%egray的最小值并輸出雙精度型
level=(fmax1-(fmax1-fmin1)/3)/255;%獲得最佳閾值
bw22=im2bw(Egray,level);%轉換圖像為二進制圖像
bw2=double(bw22);
%Step4 對得到二值圖像作開閉操作進行濾波
figure,imshow(bw2);title('圖像二值化');%得到二值圖像
grd=edge(bw2,'canny')%用canny算子識別強度圖像中的邊界
figure,imshow(grd);title('圖像邊緣提取');%輸出圖像邊緣
bg1=imclose(grd,strel('rectangle',[5,19]));%取矩形框的閉運算
figure,imshow(bg1);title('圖像閉運算[5,19]');%輸出閉運算的圖像
bg3=imopen(bg1,strel('rectangle',[5,19]));%取矩形框的開運算
figure,imshow(bg3);title('圖像開運算[5,19]');%輸出開運算的圖像
bg2=imopen(bg3,strel('rectangle',[19,1]));%取矩形框的開運算
figure,imshow(bg2);title('圖像開運算[19,1]');%輸出開運算的圖像
%Step5 對二值圖像進行區域提取，并計算區域特征參數。進行區域特征參數比較，提取車牌區域
[L,num] = bwlabel(bg2,8);%標注二進制圖像中已連接的部分
Feastats = imfeature(L,'basic');%計算圖像區域的特征尺寸
Area=[Feastats.Area];%區域面積
BoundingBox=[Feastats.BoundingBox];%[x y width height]車牌的框架大小
RGB = label2rgb(L, 'spring', 'k', 'shuffle'); %標志圖像向RGB圖像轉換
figure,imshow(RGB);title('圖像彩色標記');%輸出框架的彩色圖像
lx=0;
for l=1:num
    width=BoundingBox((l-1)*4+3);%框架寬度的計算
    hight=BoundingBox((l-1)*4+4);%框架高度的計算
    if (width>98 & width<160 & hight>25 & hight<50)%框架的寬度和高度的范圍
        lx=lx+1;
        Getok(lx)=l;
    end
end
for k= 1:lx
    l=Getok(k);    
    startcol=BoundingBox((l-1)*4+1)-2;%開始列
    startrow=BoundingBox((l-1)*4+2)-2;%開始行
    width=BoundingBox((l-1)*4+3)+8;%車牌寬
    hight=BoundingBox((l-1)*4+4)+2;%車牌高
    rato=width/hight;%計算車牌長寬比
    if rato>2 & rato<4   
        break;
    end
end
sbw1=bw2(startrow:startrow+hight,startcol:startcol+width-1); %獲取車牌二值子圖
subcol1=Sgray(startrow:startrow+hight,startcol:startcol+width-1);%獲取車牌灰度子圖
figure,subplot(2,1,1),imshow(subcol1);title('車牌灰度子圖');%輸出灰度圖像
subplot(2,1,2),imshow(sbw1);title('車牌二值子圖');%輸出車牌的二值圖
%Step6 計算車牌水平投影，并對水平投影進行峰谷分析
histcol1=sum(sbw1);      %計算垂直投影
histrow=sum(sbw1');      %計算水平投影
figure,subplot(2,1,1),bar(histcol1);title('垂直投影（含邊框）');%輸出垂直投影
subplot(2,1,2),bar(histrow);     title('水平投影（含邊框）');%輸出水平投影
figure,subplot(2,1,1),bar(histrow);     title('水平投影（含邊框）');%輸出水平投影
subplot(2,1,2),imshow(sbw1);title('車牌二值子圖');%輸出二值圖
%對水平投影進行峰谷分析
meanrow=mean(histrow);%求水平投影的平均值
minrow=min(histrow);%求水平投影的最小值
levelrow=(meanrow+minrow)/2;%求水平投影的平均值
count1=0;
l=1;
for k=1:hight
    if histrow(k)<=levelrow                             
        count1=count1+1;                                
    else 
        if count1>=1
            markrow(l)=k;%上升點
            markrow1(l)=count1;%谷寬度（下降點至下一個上升點）
            l=l+1;
        end
        count1=0;
    end
end
markrow2=diff(markrow);%峰距離（上升點至下一個上升點）
[m1,n1]=size(markrow2);
n1=n1+1;
markrow(l)=hight;
markrow1(l)=count1;
markrow2(n1)=markrow(l)-markrow(l-1);
l=0;
for k=1:n1
    markrow3(k)=markrow(k+1)-markrow1(k+1);%下降點
    markrow4(k)=markrow3(k)-markrow(k);%峰寬度（上升點至下降點）
    markrow5(k)=markrow3(k)-double(uint16(markrow4(k)/2));%峰中心位置
end 
%Step7 計算車牌旋轉角度
%(1)在上升點至下降點找第一個為1的點
[m2,n2]=size(sbw1);%sbw1的圖像大小
[m1,n1]=size(markrow4);%markrow4的大小
maxw=max(markrow4);%最大寬度為字符
if markrow4(1) ~= maxw%檢測上邊
    ysite=1;
    k1=1;
    for l=1:n2
    for k=1:markrow3(ysite)%從頂邊至第一個峰下降點掃描
        if sbw1(k,l)==1
            xdata(k1)=l;
            ydata(k1)=k;
            k1=k1+1;
            break;
        end
    end
    end
else  %檢測下邊
    ysite=n1;
    if markrow4(n1) ==0
        if markrow4(n1-1) ==maxw
           ysite= 0; %無下邊
       else
           ysite= n1-1;
       end
    end
    if ysite ~=0
        k1=1;
        for l=1:n2
            k=m2;
            while k>=markrow(ysite) %從底邊至最后一個峰的上升點掃描
                if sbw1(k,l)==1
                    xdata(k1)=l;
                    ydata(k1)=k;
                    k1=k1+1;
                    break;
                end
                k=k-1;
            end
        end
    end
end       
%(2)線性擬合，計算與x夾角
fresult = fit(xdata',ydata','poly1');   %poly1    Y = p1*x+p2
p1=fresult.p1;
angle=atan(fresult.p1)*180/pi; %弧度換為度，360/2pi,  pi=3.14
%(3)旋轉車牌圖象
subcol = imrotate(subcol1,angle,'bilinear','crop'); %旋轉車牌圖象
sbw = imrotate(sbw1,angle,'bilinear','crop');%旋轉圖像
figure,subplot(2,1,1),imshow(subcol);title('車牌灰度子圖');%輸出車牌旋轉后的灰度圖像標題顯示車牌灰度子圖
subplot(2,1,2),imshow(sbw);title('');%輸出車牌旋轉后的灰度圖像
title(['車牌旋轉角: ',num2str(angle),'度'] ,'Color','r');%顯示車牌的旋轉角度
%Step8 旋轉車牌后重新計算車牌水平投影，去掉車牌水平邊框，獲取字符高度
histcol1=sum(sbw); %計算垂直投影
histrow=sum(sbw'); %計算水平投影
figure,subplot(2,1,1),bar(histcol1);title('垂直投影（旋轉后）');
subplot(2,1,2),bar(histrow);     title('水平投影（旋轉后）');
figure,subplot(2,1,1),bar(histrow);     title('水平投影（旋轉后）');
subplot(2,1,2),imshow(sbw);title('車牌二值子圖（旋轉后）');
%去水平（上下）邊框,獲取字符高度
maxhight=max(markrow2);
findc=find(markrow2==maxhight);
rowtop=markrow(findc);
rowbot=markrow(findc+1)-markrow1(findc+1);
sbw2=sbw(rowtop:rowbot,:);  %子圖為(rowbot-rowtop+1)行
maxhight=rowbot-rowtop+1;   %字符高度(rowbot-rowtop+1)
%Step9 計算車牌垂直投影，去掉車牌垂直邊框，獲取車牌及字符平均寬度
histcol=sum(sbw2);  %計算垂直投影
figure,subplot(2,1,1),bar(histcol);title('垂直投影（去水平邊框后）');%輸出車牌的垂直投影圖像
subplot(2,1,2),imshow(sbw2); %輸出垂直投影圖像
title(['車牌字符高度： ',int2str(maxhight)],'Color','r');%輸出車牌字符高度
%對垂直投影進行峰谷分析
meancol=mean(histcol);%求垂直投影的平均值
mincol=min(histcol);%求垂直投影的平均值
levelcol=(meancol+mincol)/4;%求垂直投影的1/4
count1=0;
l=1;
for k=1:width
    if histcol(k)<=levelcol 
        count1=count1+1;
    else 
        if count1>=1
            markcol(l)=k; %字符上升點
            markcol1(l)=count1; %谷寬度（下降點至下一個上升點）
            l=l+1;
        end
        count1=0;
    end
end
markcol2=diff(markcol);%字符距離（上升點至下一個上升點）
[m1,n1]=size(markcol2);
n1=n1+1;
markcol(l)=width;
markcol1(l)=count1;
markcol2(n1)=markcol(l)-markcol(l-1);
%Step10 計算車牌上每個字符中心位置，計算最大字符寬度maxwidth
l=0;
for k=1:n1
    markcol3(k)=markcol(k+1)-markcol1(k+1);%字符下降點
    markcol4(k)=markcol3(k)-markcol(k); %字符寬度（上升點至下降點）
    markcol5(k)=markcol3(k)-double(uint16(markcol4(k)/2));%字符中心位置
end 
markcol6=diff(markcol5); %字符中心距離（字符中心點至下一個字符中心點）
maxs=max(markcol6); %查找最大值，即為第二字符與第三字符中心距離
findmax=find(markcol6==maxs);
markcol6(findmax)=0;
maxwidth=max(markcol6);%查找最大值，即為最大字符寬度
%Step11 提取分割字符,并變換為22行*14列標準子圖
l=1;
[m2,n2]=size(subcol);
figure;
for k=findmax-1:findmax+5
        cleft=markcol5(k)-maxwidth/2;
        cright=markcol5(k)+maxwidth/2-2;
        if cleft<1
            cleft=1;
            cright=maxwidth;
        end
        if cright>n2
            cright=n2;
            cleft=n2-maxwidth;
        end
        SegGray=subcol(rowtop:rowbot,cleft:cright);
        SegBw1=sbw(rowtop:rowbot,cleft:cright);
        SegBw2 = imresize(SegBw1,[22 14]);              %變換為32行*16列標準子圖      
        subplot(2,n1,l),imshow(SegGray);
        if l==7
            title(['車牌字符寬度： ',int2str(maxwidth)],'Color','r');
        end
        subplot(2,n1,n1+l),imshow(SegBw2);        
        fname=strcat('F:\MATLAB\work\sam\image',int2str(k),'.jpg');%保存子圖備選入樣本庫,并建立樣本庫
        imwrite(SegBw2,fname,'jpg') 
        l=l+1;
end
%Step12 將計算計算獲取的字符圖像與樣本庫進行匹配，自動識別出字符代碼。
%在進行匹配前需將存入的子圖名按以下編碼規則更改：
%Sam0 ~Sam9   0~9 數字  1~10
%SamA~SamZ  A~Z 字母    11~36
%Sam粵~Sam滬  省別 漢字 37~45
liccode=char(['0':'9' 'A':'Z' '粵桂海云貴川京津滬']); %建立自動識別字符代碼表  
SubBw2=zeros(22,14);
l=1;
[m2,n2]=size(subcol);
for k=findmax-1:findmax+5
       cleft=markcol5(k)-maxwidth/2;
        cright=markcol5(k)+maxwidth/2-2;
        if cleft<1
            cleft=1;
            cright=maxwidth;
        end
        if cright>n2
            cright=n2;
            cleft=n2-maxwidth;
        end
        SegBw1=sbw(rowtop:rowbot,cleft:cright);
        SegBw2 = imresize(SegBw1,[22 14]);%變換為22行*14列標準子圖      
        if l==1                 %第一位漢字識別
            kmin=37;
            kmax=45;
        elseif l==2             %第二位 A~Z 字母識別
            kmin=11;
            kmax=36;
        elseif l>=3 & l<=5      %第三、四位 0~9  A~Z字母和數字識別
            kmin=1;
            kmax=36;
        else                    %第五～七位 0~9 數字識別
            kmin=1;
            kmax=10;
        end
        for k2=kmin:kmax
            fname=strcat('F:\MATLAB\work\sam\Sam',liccode(k2),'.jpg');
            SamBw2 = imread(fname);
            %SamBw2 =double(SamBw2);
            %SubBw2 = SamBw2-SegBw2;
            for  i=1:32
                for j=1:16
                    SubBw2(i,j)=SegBw2(i,j)-SamBw2(i,j);
                end
            end
            Dmax=0;
            for k1=1:32
                for l1=1:16
                    if ( SubBw2(k1,l1) > 0 | SubBw2(k1,l1) <0 )
                        Dmax=Dmax+1;
                    end
                end
            end
            Error(k2)=Dmax;
        end
        Error1=Error(kmin:kmax);
        MinError=min(Error1);
        findc=find(Error1==MinError);
        RegCode(l*2-1)=liccode(findc(1)+kmin-1);
        RegCode(l*2)=' ';
        l=l+1;
end
title (['識別車牌號碼:', RegCode],'Color','r');