% FREQEST - Estimate fingerprint ridge frequency within image block

%

% Function to estimate the fingerprint ridge frequency within a small block

% of a fingerprint image.  This function is used by RIDGEFREQ

%

% Usage:

%  freqim =  freqest(im, orientim, windsze, minWaveLength, maxWaveLength)

%

% Arguments:

%         im       - Image block to be processed.

%         orientim - Ridge orientation image of image block.

%         windsze  - Window length used to identify peaks. This should be

%                    an odd integer, say 3 or 5.

%         minWaveLength,  maxWaveLength - Minimum and maximum ridge

%                     wavelengths, in pixels, considered acceptable.

% 

% Returns:

%         freqim    - An image block the same size as im with all values

%                     set to the estimated ridge spatial frequency.  If a

%                     ridge frequency cannot be found, or cannot be found

%                     within the limits set by min and max Wavlength

%                     freqim is set to zeros.

%

% Suggested parameters for a 500dpi fingerprint image

%   freqim = freqest(im,orientim, 5, 5, 15);

%

% See also:  RIDGEFREQ, RIDGEORIENT, RIDGESEGMENT

%

% Note I am not entirely satisfied with the output of this function.

% Peter Kovesi 

% School of Computer Science & Software Engineering

% The University of Western Australia

% pk at csse uwa edu au

% http://www.csse.uwa.edu.au/~pk

%

% January 2005

function freqim =  freqest(im, orientim, windsze, minWaveLength, maxWaveLength)

    debug = 0;

    [rows,cols] = size(im);

    % Find mean orientation within the block. This is done by averaging the

    % sines and cosines of the doubled angles before reconstructing the

    % angle again.  This avoids wraparound problems at the origin.

    orientim = 2*orientim(:);    

    cosorient = mean(cos(orientim));

    sinorient = mean(sin(orientim));    

    orient = atan2(sinorient,cosorient)/2;

    % Rotate the image block so that the ridges are vertical

    rotim = imrotate(im,orient/pi*180+90,'nearest', 'crop');

    % Now crop the image so that the rotated image does not contain any

    % invalid regions.  This prevents the projection down the columns

    % from being mucked up.

    cropsze = fix(rows/sqrt(2)); offset = fix((rows-cropsze)/2);

    rotim = rotim(offset:offset+cropsze, offset:offset+cropsze);

    % Sum down the columns to get a projection of the grey values down

    % the ridges.

    proj = sum(rotim);

    % Find peaks in projected grey values by performing a greyscale

    % dilation and then finding where the dilation equals the original

    % values. 

    dilation = ordfilt2(proj, windsze, ones(1,windsze));

    maxpts = (dilation == proj) & (proj > mean(proj));

    maxind = find(maxpts);

    % Determine the spatial frequency of the ridges by divinding the

    % distance between the 1st and last peaks by the (No of peaks-1). If no

    % peaks are detected, or the wavelength is outside the allowed bounds,

    % the frequency image is set to 0

    if length(maxind) < 2

    freqim = zeros(size(im));

    else

    NoOfPeaks = length(maxind);

    waveLength = (maxind(end)-maxind(1))/(NoOfPeaks-1);

    if waveLength > minWaveLength & waveLength < maxWaveLength

        freqim = 1/waveLength * ones(size(im));

    else

        freqim = zeros(size(im));

    end

    end

    if debug

    show(im,1)

    show(rotim,2);

    figure(3),    plot(proj), hold on

    meanproj = mean(proj)

    if length(maxind) < 2

        fprintf('No peaks found\n');

    else

        plot(maxind,dilation(maxind),'r*'), hold off

        waveLength = (maxind(end)-maxind(1))/(NoOfPeaks-1);

    end

    end