%根據(jù)連續(xù)分布方向圖求塊方向圖 利用8*8的窗口
clear;close all;
fp=imread('105_2.tif');  %讀入指紋圖象
imshow(fp);title('原始指紋圖象')               %顯示指紋圖象 
M0=80;  
Var0=200;
[len,wid]=size(fp);
fpn=zeros(len,wid);       %為規(guī)格化后的圖象設置存儲空間
M=mean2(fp);              %求原始指紋圖象灰度均值
Var=std2(fp)^2;             %求原始指紋圖象灰度方差
%規(guī)格化
fp=double(fp);   %將unit8型轉換成double，便于sqrt的計算
for i=1:len
    for j=1:wid
        if fp(i,j)>M
            fpn(i,j)=M0+sqrt( Var0*(fp(i,j)-M)^2/Var );  
        else
            fpn(i,j)=M0-sqrt( Var0*(fp(i,j)-M)^2/Var );
        end
    end
end
fp=uint8(fpn);
figure,imshow(fp);title('規(guī)格化后的指紋圖象');  %顯示規(guī)格化后的指紋圖象
block_map=zeros(len,wid);
fp=double(fp);
%利用9*9的掩模，邊緣的四個像素不能處理
for i=5:len-4
    for j=5:wid-4                            
         %求每個像素在8個方向上的灰度平均值；
        aver_gray_d1=sum( [fp(i,j),fp(i,j-4),fp(i,j-2),fp(i,j+2),fp(i,j+4)] )/5;                   
        aver_gray_d2=sum([fp(i,j),fp(i-2,j-4),fp(i-1,j-2),fp(i+1,j+2),fp(i+2,j+4)])/5;
        aver_gray_d3=sum([fp(i,j),fp(i-4,j-4),fp(i-2,j-2),fp(i+2,j+2),fp(i+4,j+4)])/5;
        aver_gray_d4=sum([fp(i,j),fp(i-4,j-2),fp(i-2,j-1),fp(i+2,j+1),fp(i+4,j+2)])/5;
        aver_gray_d5=sum([fp(i,j),fp(i-4,j),fp(i-2,j),fp(i+2,j),fp(i+4,j)])/5;
        aver_gray_d6=sum([fp(i,j),fp(i+4,j-2),fp(i+2,j-1),fp(i-2,j+1),fp(i-4,j+2)])/5;
        aver_gray_d7=sum([fp(i,j),fp(i+4,j-4),fp(i+2,j-2),fp(i-2,j+2),fp(i-4,j+4)])/5;
        aver_gray_d8=sum([fp(i,j),fp(i+2,j-4),fp(i+1,j-2),fp(i-1,j+2),fp(i-2,j+4)])/5;
         %求每個像素在8個方向上的方差；
        std_d1=sqrt( sum( [ (aver_gray_d1-fp(i,j-4))^2,(aver_gray_d1-fp(i,j-2))^2,(aver_gray_d1-fp(i,j+2))^2,(aver_gray_d1-fp(i,j+4))^2 ] )/4);
        std_d2=sqrt( sum( [(aver_gray_d2-fp(i-2,j-4))^2,(aver_gray_d2-fp(i-1,j-2))^2,(aver_gray_d2-fp(i+1,j+2))^2,(aver_gray_d2-fp(i+2,j+4))^2] )/4);
        std_d3=sqrt( sum( [(aver_gray_d3-fp(i-4,j-4))^2,(aver_gray_d3-fp(i-2,j-2))^2,(aver_gray_d3-fp(i+2,j+2))^2,(aver_gray_d3-fp(i+4,j+4))^2] )/4);
        std_d4=sqrt( sum( [(aver_gray_d4-fp(i-4,j-2))^2,(aver_gray_d4-fp(i-2,j-1))^2,(aver_gray_d4-fp(i+2,j+1))^2,(aver_gray_d4-fp(i+4,j+2))^2] )/4);
        std_d5=sqrt( sum( [(aver_gray_d5-fp(i-4,j))^2,(aver_gray_d5-fp(i-2,j))^2,(aver_gray_d5-fp(i+2,j))^2,(aver_gray_d5-fp(i+4,j))^2] )/4);
        std_d6=sqrt( sum( [(aver_gray_d6-fp(i+4,j-2))^2,(aver_gray_d6-fp(i+2,j-1))^2,(aver_gray_d6-fp(i-2,j+1))^2,(aver_gray_d6-fp(i-4,j+2))^2] )/4);
        std_d7=sqrt( sum( [(aver_gray_d7-fp(i+4,j-4))^2,(aver_gray_d7-fp(i+2,j-2))^2,(aver_gray_d7-fp(i-2,j+2))^2,(aver_gray_d7-fp(i-4,j+4))^2] )/4);
        std_d8=sqrt( sum( [(aver_gray_d8-fp(i+2,j-4))^2,(aver_gray_d8-fp(i+1,j-2))^2,(aver_gray_d8-fp(i-1,j+2))^2,(aver_gray_d8-fp(i-2,j+4))^2] )/4);
        std_d=...
            [std_d1,std_d2,std_d3,std_d4,std_d5,std_d6,std_d7,std_d8];
        [minstd,minstd_index]=min(std_d);
        block_map(i,j)=minstd_index;
    end
end

%下面顯示點方向圖，利用索引圖象，imshow(bit_dmap,map)
%去除邊緣的像素,得到點方向圖的數(shù)據(jù)矩陣bit_dmap
bit_dmap=imcrop(block_map,[5 5 247 247]);        %去除邊緣的4個像素  
%定義map函數(shù)
map=[1 0 0;          %方向1為紅色
     0 1 0;          %方向2為綠色
     0 0 1;          %方向3為藍色
     1 1 0;          %方向4為紅+綠=黃色
     1 0 1;          %方向5為紅+藍=洋紅色
     0 1 1;          %方向6為綠+藍=青色
     1 1 1;          %方向7為白色
     0 0 0];         %方向8為黑色
figure,imshow(bit_dmap,map);colorbar;
title('改進方法得到的點方向圖');
text(300,80,{'方向8為22.5度' ,'方向7為45度','方向6為67.5度','方向5為-90度','方向4為-67.5度','方向3為-45度','方向2為-22.5度','方向1為0度'});

%下面通過點方向圖求解塊連續(xù)分布方向圖
%原理：采用平滑濾波技術，對點方向圖進行平滑處理。
cd_dmap=zeros(len,wid);
for  i=5:len-4
  for       j=5:wid-4                    
         blocknum=zeros(8,1);
         for m=-4:4                %取9×9的模板
            for n=-4:4
                 
                 switch  block_map(i-m,j-n)  %在模板內計算方向直方圖
                       case 1   
                             blocknum(1)=blocknum(1)+1;           %方向1的像素數(shù)目
                       case 2   
                             blocknum(2)=blocknum(2)+1;
                       case 3   
                             blocknum(3)=blocknum(3)+1;
                       case 4   
                             blocknum(4)=blocknum(4)+1;
                       case 5   
                             blocknum(5)=blocknum(5)+1;
                       case 6   
                             blocknum(6)=blocknum(6)+1;
                       case 7   
                             blocknum(7)=blocknum(7)+1;
                       case 8   
                             blocknum(8)=blocknum(8)+1;
                      end
           end                                                         %模板的最大方向確定
       end
      [cd_maxd,cd_max_index]=max(blocknum);     %最大的方向賦給block_maxd_index，該方向的像素數(shù)目為cd_max_index
      cd_dmap(i,j)=cd_max_index;
  end
end
figure,imshow(imcrop(cd_dmap,[5 5 247 247]),map);colorbar;
title('連續(xù)分布方向圖');
text(300,80,{'方向8為22.5度' ,'方向7為45度','方向6為67.5度','方向5為-90度','方向4為-67.5度','方向3為-45度','方向2為-22.5度','方向1為0度'});



%以下計算塊方向圖
%先將圖象分塊，通用做法分成32*32塊,每塊為8*8
block_dmap=ones(256,256);                   %給塊方向圖定義存儲空間，初值為1，全白色
for  i=1:32
  for       j=1:32                   %分成32*32個塊
       x=cd_dmap([1+(i-1)*8:8+(i-1)*8],[1+(j-1)*8:8+(j-1)*8]);  %未做平滑       
       y=block_dmap([1+(i-1)*8:8+(i-1)*8],[1+(j-1)*8:8+(j-1)*8]); %未做平滑 
       y_size=size(y);
       blocknum=zeros(8,1);
       for m=1:8                %每個小塊由8*8像素組成
           for n=1:8
                 
                 switch  x(m,n)  %在每小塊內計算方向直方圖
                       case 1   
                             blocknum(1)=blocknum(1)+1;           %方向1的像素數(shù)目
                       case 2   
                             blocknum(2)=blocknum(2)+1;
                       case 3   
                             blocknum(3)=blocknum(3)+1;
                       case 4   
                             blocknum(4)=blocknum(4)+1;
                       case 5   
                             blocknum(5)=blocknum(5)+1;
                       case 6   
                             blocknum(6)=blocknum(6)+1;
                       case 7   
                             blocknum(7)=blocknum(7)+1;
                       case 8   
                             blocknum(8)=blocknum(8)+1;

                 end
           end                                                         %每個小塊的方向確定
       end
                [block_maxd,block_maxd_index]=max(blocknum);     %最大的方向賦給block_maxd_index，該方向的像素數(shù)目為block_maxd
                %根據(jù)block_maxd_index的值，在該小塊內繪制該方向
                switch block_maxd_index
                case 1
                  y(4,[2:7])=0;                 %方向1，即方向為0度時，令第8行第2-15列的元素為0（黑色）
                case 2
                  idx=sub2ind(y_size,[  3 4 4 5 5 6  ],[2:7]);
                  y(idx)=0;                 
                case 3
                  idx=sub2ind(y_size,[2:7],[2:7]);
                  y(idx)=0;
                case 4
                  idx=sub2ind(y_size,[2:7],[  3 4 4 5 5 6  ]);
                  y(idx)=0;
                case 5
                  y([2:7],4)=0;
                case 6
                  idx=sub2ind(y_size,[7:-1:2],[   3 4 4 5 5 6 ]);
                  y(idx)=0;
                case 7
                  idx=sub2ind(y_size,[7:-1:2],[7:-1:2]);
                  y(idx)=0;                     
                case 8
                  idx=sub2ind(y_size,[  6 5 5 4 4 3  ],[2:7]);
                  y(idx)=0;
                end
                block_dmap([1+(i-1)*8:8+(i-1)*8],[1+(j-1)*8:8+(j-1)*8])=y;
          end
end
figure;imshow(block_dmap);title('由指紋連續(xù)分布方向圖求得的塊方向圖');
fp=double(fp);
rever_block_dmap=1-block_dmap;               %rever_block_dmap內的元素非0即1
rever_block_dmap=50*rever_block_dmap;       %將rever_block_dmap內的灰度值乘上50，以便顯示疊加的圖形
addblock_dmap=imadd(fp,rever_block_dmap,'double');
figure;imshow(addblock_dmap,[]); title('原指紋圖象和塊方向圖的疊加圖形');     %顯示原指紋圖象和塊方向圖疊加的圖形