% Bc. Michal Kubiš
% Fakulta elektrotechniky a komunikačních technologií na VUT v Brně
% Vedoucí: Ing. Petr Fedra
% 
% 21. 5. 2010
% založeno na algoritmu Petera Kovesiho

%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
function [oimg,fimg,bwimg,eimg,enhimg] =  fft_enhance_cubs(img)
    global NFFT;
    
    if evalin('base','exist(''OVRLP'')') == 0     
        assignin('base', 'NFFT', 32)               
        assignin('base', 'BLKSZ', 12)                  
        assignin('base', 'OVRLP', 6)
        assignin('base', 'ALPHA', 0.4)
        assignin('base', 'RMIN', 5)
        assignin('base', 'RMAX', 20)
    end

    NFFT = evalin('base', 'NFFT');        
    BLKSZ = evalin('base', 'BLKSZ');
    OVRLP = evalin('base', 'OVRLP');
    ALPHA = evalin('base', 'ALPHA');
    RMIN = evalin('base', 'RMIN');
    RMAX = evalin('base', 'RMAX');


%     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