%分組法求點方向圖
clear;close all;
fp=imread('finger.tif');  %讀入指紋圖象
imshow(fp);title('原始指紋圖象')               %顯示指紋圖象  
%以下是對原始指紋圖象進行規格化，規格化后的灰度均值M0設為80，方差Var0為200
M0=80;  
Var0=200;
[len,wid]=size(fp);
fpn=zeros(len,wid);       %為規格化后的圖象設置存儲空間
M=mean2(fp);              %求原始指紋圖象灰度均值
Var=std2(fp)^2;             %求原始指紋圖象灰度方差
%規格化
fp=double(fp);   %將uint8型轉換成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
fpn=uint8(fpn);
figure,imshow(fpn);title('規格化后的指紋圖象');  %顯示規格化后的指紋圖象
aver_d=zeros(len,wid);
%利用9*9的掩模，邊緣的四個像素不能處理
for i=5:len-4
    for j=5:wid-4                            
        aver_gray_d1=sum([fpn(i,j-4),fpn(i,j-2),fpn(i,j+2),fpn(i,j+4)])/4;                    %求8個方向上的灰度平均值；
        aver_gray_d2=sum([fpn(i-2,j-4),fpn(i-1,j-2),fpn(i+1,j+2),fpn(i+2,j+4)])/4;
        aver_gray_d3=sum([fpn(i-4,j-4),fpn(i-2,j-2),fpn(i+2,j+2),fpn(i+4,j+4)])/4;
        aver_gray_d4=sum([fpn(i-4,j-2),fpn(i-2,j-1),fpn(i+2,j+1),fpn(i+4,j+2)])/4;
        aver_gray_d5=sum([fpn(i-4,j),fpn(i-2,j),fpn(i+2,j),fpn(i+4,j)])/4;
        aver_gray_d6=sum([fpn(i+4,j-2),fpn(i+2,j-1),fpn(i-2,j+1),fpn(i-4,j+2)])/4;
        aver_gray_d7=sum([fpn(i+4,j-4),fpn(i+2,j-2),fpn(i-2,j+2),fpn(i-4,j+4)])/4;
        aver_gray_d8=sum([fpn(i+2,j-4),fpn(i+1,j-2),fpn(i-1,j+2),fpn(i-2,j+4)])/4;
        d_diff1=abs(aver_gray_d1-aver_gray_d5);                                    %對兩兩垂直方向的灰度平均值差求絕對值；
        d_diff2=abs(aver_gray_d2-aver_gray_d6);
        d_diff3=abs(aver_gray_d3-aver_gray_d7);
        d_diff4=abs(aver_gray_d4-aver_gray_d8);
        [max_d,max_d_index]=max([d_diff1,d_diff2,d_diff3,d_diff4]);           %取差值絕對值最大的兩個方向作為可能的脊線方向；
        %方向max_d_index和max_d_index+4為fp(i,j)可能的脊線方向
                                                                                           
        aver_gray_d=...
          [aver_gray_d1,aver_gray_d2,aver_gray_d3,aver_gray_d4,aver_gray_d5,aver_gray_d6,aver_gray_d7,aver_gray_d8];
        %取max_d_index和max_d_index+4兩個方向中灰度平均值與fp(i,j)比較接近的方向，作為點（i,j）的方向；
            if abs(fpn(i,j)-aver_gray_d(max_d_index))<abs(fpn(i,j)-aver_gray_d(max_d_index+4))
                  aver_d(i,j)=max_d_index;
            else
                  aver_d(i,j)=max_d_index+4;
            end
    end
end

%下面顯示點方向圖，利用索引圖象，imshow(bit_dmap,map)
%去除邊緣的像素,得到點方向圖的數據矩陣bit_dmap
bit_dmap=imcrop(aver_d,[5 5 247 247]);        %去除邊緣的4個像素  
%定義d_map函數
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,{'利用9*9模板，計算4組','兩兩垂直方向的','灰度平均值','所得到的點方向圖',...
   '其中','方向8為22.5度' ,'方向7為45度','方向6為67.5度','方向5為-90度','方向4為-67.5度','方向3為-45度','方向2為-22.5度','方向1為0度'});