% 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