fft_enhance_cubs.m

指纹增强算法, 研究指纹识别的可以参考一下. matlab.

%--------------------------------------------------------------------------
%fft_enhance_cubs
%enhances the fingerprint image
%syntax:
%[oimg,fimg,bwimg,eimg,enhimg] =  fft_enhance_cubs(img)
%oimg -  [OUT] block orientation image(can be viewed using
%        view_orientation_image.m)
%fimg  - [OUT] block frequency image(indicates ridge spacing)
%bwimg - [OUT] shows angular bandwidth image(filter bandwidth adapts near the
%        singular points)
%eimg  - [OUT] energy image. Indicates the 'ridgeness' of a block (can be 
%        used for fingerprint segmentation)
%enhimg- [OUT] enhanced image
%img   - [IN]  input fingerprint image (HAS to be of DOUBLE type)
%Contact:
%   [e-mail] [e-mail]
%   http://www.sharat.org
%Reference:
%1. S. Chikkerur, C.Wu and V. Govindaraju, "Systematic approach for feature
%   extraction in Fingerprint Images", ICBA 2004
%2. S. Chikkerur and V. Govindaraju, "Fingerprint Image Enhancement using 
%   STFT Analysis", International Workshop on Pattern Recognition for Crime 
%   Prevention, Security and Surveillance, ICAPR 2005
%3. S. Chikkeur, "Online Fingerprint Verification", M. S. Thesis,
%   University at Buffalo, 2005
%4. T. Jea and V. Govindaraju, "A Minutia-Based Partial Fingerprint Recognition System", 
%   to appear in Pattern Recognition 2005
%5. S. Chikkerur, "K-plet and CBFS: A Graph based Fingerprint
%   Representation and Matching Algorithm", submitted, ICB 2006
% See also: cubs_visualize_template
%--------------------------------------------------------------------------
function [oimg,fimg,bwimg,eimg,enhimg] =  fft_enhance_cubs(img)
    global NFFT;
    NFFT            =   32;     %size of FFT
    BLKSZ           =   12;      %size of the block
    OVRLP           =   6;      %size of overlap
    ALPHA           =   0.4;    %root filtering
    RMIN            =   5;      %min allowable ridge spacing
    RMAX            =   20;     %maximum allowable ridge spacing
    do_prefiltering =   1;
    [nHt,nWt]       =   size(img);  
    img             =   im2double(img);    %convert to DOUBLE
    
    nBlkHt          =   floor((nHt-2*OVRLP)/BLKSZ);
    nBlkWt          =   floor((nWt-2*OVRLP)/BLKSZ);
    fftSrc          =   zeros(nBlkHt*nBlkWt,NFFT*NFFT); %stores FFT
    nWndSz          =   BLKSZ+2*OVRLP; %size of analysis window. 
    warning off MATLAB:divideByZero
    %-------------------------
    %allocate outputs
    %-------------------------
    oimg        =   zeros(nBlkHt,nBlkWt);
    fimg        =   zeros(nBlkHt,nBlkWt);
    bwimg       =   zeros(nBlkHt,nBlkWt);
    eimg        =   zeros(nBlkHt,nBlkWt);
    enhimg      =   zeros(nHt,nWt);
    
    %-------------------------
    %precomputations
    %-------------------------
    [x,y]       =   meshgrid(0:nWndSz-1,0:nWndSz-1);
    dMult       =   (-1).^(x+y); %used to center the FFT
    [x,y]       =   meshgrid(-NFFT/2:NFFT/2-1,-NFFT/2:NFFT/2-1);
    r           =   sqrt(x.^2+y.^2)+eps;
    th          =   atan2(y,x);
    th(th<0)    =   th(th<0)+pi;
    w           =   raised_cosine_window(BLKSZ,OVRLP); %spectral window
    
    %-------------------------
    %FFT Analysis
    %-------------------------
    for i = 0:nBlkHt-1
        nRow = i*BLKSZ+OVRLP+1;  
        for j = 0:nBlkWt-1
            nCol = j*BLKSZ+OVRLP+1;
            %extract local block
            blk     =   img(nRow-OVRLP:nRow+BLKSZ+OVRLP-1,nCol-OVRLP:nCol+BLKSZ+OVRLP-1);
            %remove dc
            dAvg    =   sum(sum(blk))/(nWndSz*nWndSz);
            blk     =   blk-dAvg;   %remove DC content
            blk     =   blk.*w;     %multiply by spectral window
            %--------------------------
            %do pre filtering
            %--------------------------
            blkfft  =   fft2(blk.*dMult,NFFT,NFFT);
            if(do_prefiltering)
                dEnergy =   abs(blkfft).^2;
                blkfft  =   blkfft.*sqrt(dEnergy);      %root filtering(for diffusion)
            end;
            fftSrc(nBlkWt*i+j+1,:) = transpose(blkfft(:));
            dEnergy =   abs(blkfft).^2;             %----REDUCE THIS COMPUTATION----
            %--------------------------
            %compute statistics
            %--------------------------
            dTotal          =   sum(sum(dEnergy));%/(NFFT*NFFT);
            fimg(i+1,j+1)   =   NFFT/(compute_mean_frequency(dEnergy,r)+eps); %ridge separation
            oimg(i+1,j+1)   =   compute_mean_angle(dEnergy,th);         %ridge angle
            eimg(i+1,j+1)   =   log(dTotal+eps);                        %used for segmentation
        end;%for j
    end;%for i

    %-------------------------
    %precomputations
    %-------------------------
    [x,y]       =   meshgrid(-NFFT/2:NFFT/2-1,-NFFT/2:NFFT/2-1);
    dMult       =   (-1).^(x+y); %used to center the FFT

    %-------------------------
    %process the resulting maps
    %-------------------------
    for i = 1:3
        oimg = smoothen_orientation_image(oimg);            %smoothen orientation image
    end;
    fimg    =   smoothen_frequency_image(fimg,RMIN,RMAX,5); %diffuse frequency image
    cimg    =   compute_coherence(oimg);                    %coherence image for bandwidth
    bwimg   =   get_angular_bw_image(cimg);                 %QUANTIZED bandwidth image
    %-------------------------
    %FFT reconstruction
    %-------------------------
    for i = 0:nBlkHt-1
        for j = 0:nBlkWt-1
            nRow = i*BLKSZ+OVRLP+1;            
            nCol = j*BLKSZ+OVRLP+1;
            %--------------------------
            %apply the filters
            %--------------------------
            blkfft  =   reshape(transpose(fftSrc(nBlkWt*i+j+1,:)),NFFT,NFFT);
            %--------------------------
            %reconstruction
            %--------------------------
            af      =   get_angular_filter(oimg(i+1,j+1),bwimg(i+1,j+1));
            rf      =   get_radial_filter(fimg(i+1,j+1));
            blkfft  =   blkfft.*(af).*(rf); 
            blk     =   real(ifft2(blkfft).*dMult);
            enhimg(nRow:nRow+BLKSZ-1,nCol:nCol+BLKSZ-1)=blk(OVRLP+1:OVRLP+BLKSZ,OVRLP+1:OVRLP+BLKSZ);
        end;%for j
    end;%for i
    %end block processing
    %--------------------------
    %contrast enhancement
    %--------------------------
    enhimg =sqrt(abs(enhimg)).*sign(enhimg);
    enhimg =imscale(enhimg);
    enhimg =im2uint8(enhimg);
    
    %--------------------------
    %clean up the image
    %--------------------------
    emsk            = segment_print(enhimg,0);
    enhimg(emsk==0) = 128;
%end function fft_enhance_cubs

%-----------------------------------
%raised_cosine
%returns 1D spectral window
%syntax:
%y = raised_cosine(nBlkSz,nOvrlp)
%y      - [OUT] 1D raised cosine function
%nBlkSz - [IN]  the window is constant here
%nOvrlp - [IN]  the window has transition here
%-----------------------------------
function y = raised_cosine(nBlkSz,nOvrlp)
    nWndSz  =   (nBlkSz+2*nOvrlp);
    x       =   abs(-nWndSz/2:nWndSz/2-1);
    y       =   0.5*(cos(pi*(x-nBlkSz/2)/nOvrlp)+1);
    y(abs(x)<nBlkSz/2)=1;
%end function raised_cosine

%-----------------------------------
%raised_cosine_window
%returns 2D spectral window
%syntax:
%w = raised_cosine_window(blksz,ovrlp)
%w      - [OUT] 1D raised cosine function
%nBlkSz - [IN]  the window is constant here
%nOvrlp - [IN]  the window has transition here
%-----------------------------------
function w = raised_cosine_window(blksz,ovrlp)
    y = raised_cosine(blksz,ovrlp);
    w = y(:)*y(:)';
%end function raised_cosine_window

%---------------------------------------------------------------------
%get_angular_filter
%generates an angular filter centered around 'th' and with bandwidth 'bw'
%the filters in angf_xx are precomputed using angular_filter_bank.m
%syntax:
%r = get_angular_filter(t0,bw)
%r - [OUT] angular filter of size NFFTxNFFT
%t0- mean angle (obtained from orientation image)
%bw- angular bandwidth(obtained from bandwidth image)
%angf_xx - precomputed filters (using angular_filter_bank.m)
%-----------------------------------------------------------------------
function rmsk = get_angular_filter(t0,BW)
     global NFFT;
     [x,y]   =   meshgrid(-NFFT/2:NFFT/2-1,-NFFT/2:NFFT/2-1);
     r       =   sqrt(x.^2+y.^2);
     th      =   atan2(y,x);
     th(th<0)=   th(th<0)+2*pi;  %unsigned
     t1      =   mod(t0+pi,2*pi);
     %-----------------
     %first lobe
     %-----------------
     d          = angular_distance(th,t0);
     msk        = 1+cos(d*pi/BW); 
     msk(d>BW)  = 0;
     rmsk       = msk;                              %save first lobe

     %-----------------
     %second lobe
     %-----------------
     d          = angular_distance(th,t1);
     msk        = 1+cos(d*pi/BW); 
     msk(d>BW)  = 0;
     rmsk       = (rmsk+msk);
%end function



%---------------------------------------------------------------------
%get_radial_filter
%generates an radial filter
%syntax:
%r = get_radial_filter(r0)
%r - [OUT] angular filter of size NFFTxNFFT
%r0- center frequency
%-----------------------------------------------------------------------
function rmsk = get_radial_filter(r0)
     global NFFT;
     N          =   4;
     r0         =   NFFT/r0;
     BW         =   r0*1.75-r0/1.75;
     [x,y]      =   meshgrid(-NFFT/2:NFFT/2-1,-NFFT/2:NFFT/2-1);
     r          =   sqrt(x.^2+y.^2);
     num        =   (r*BW).^(2*N);
     den        =   (r*BW).^(2*N)+((r.^2-r0.^2)).^(2*N);
     rmsk       =   sqrt(num./den);
%end function get_radial_filter

%-----------------------------------------------------------
%get_angular_bw_image
%the bandwidth allocation is currently based on heuristics
%(domain knowledge :)). 
%syntax:
%bwimg = get_angular_bw_image(c)
%-----------------------------------------------------------
function bwimg = get_angular_bw_image(c)
    bwimg   =   zeros(size(c));
    bwimg(:,:)    = pi/2;                       %med bw
    bwimg(c<=0.7) = pi;                         %high bw
    bwimg(c>=0.9) = pi/6;                       %low bw
%end function get_angular_bw


%-----------------------------------------------------------
%get_angular_bw_image
%the bandwidth allocation is currently based on heuristics
%(domain knowledge :)). 
%syntax:
%bwimg = get_angular_bw_image(c)
%-----------------------------------------------------------
function mth = compute_mean_angle(dEnergy,th)
    global NFFT;
    sth         =   sin(2*th);
    cth         =   cos(2*th);
    num         =   sum(sum(dEnergy.*sth));
    den         =   sum(sum(dEnergy.*cth));
    mth         =   0.5*atan2(num,den);
    if(mth <0)
        mth = mth+pi;
    end;
%end function compute_mean_angle

%-----------------------------------------------------------
%get_angular_bw_image
%the bandwidth allocation is currently based on heuristics
%(domain knowledge :)). 
%syntax:
%bwimg = get_angular_bw_image(c)
%-----------------------------------------------------------
function mr = compute_mean_frequency(dEnergy,r)
    global NFFT;
    num         =   sum(sum(dEnergy.*r));
    den         =   sum(sum(dEnergy));
    mr          =   num/(den+eps);
%end function compute_mean_angle


%-----------------------------------------
%angular_distance
%computes angular distance-acute angle 
%-----------------------------------------
function d = angular_distance(th,t0)
    d = abs(th-t0);
    d = min(d,2*pi-d);
%end function angular_distance