% 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


    
function freqim =  freqest(im, orientim, windsze, minWaveLength, maxWaveLength)
    
    debug = 0;
    
    [rows,cols] = size(im);
    
    orientim = 2*orientim(:);    
    cosorient = mean(cos(orientim));
    sinorient = mean(sin(orientim));    
    orient = atan2(sinorient,cosorient)/2;

    rotim = imrotate(im,orient/pi*180+90,'nearest', 'crop');
    
    cropsze = fix(rows/sqrt(2)); offset = fix((rows-cropsze)/2);
    rotim = rotim(offset:offset+cropsze, offset:offset+cropsze);

    proj = sum(rotim);
    
    dilation = ordfilt2(proj, windsze, ones(1,windsze));
    maxpts = (dilation == proj) & (proj > mean(proj));
    maxind = find(maxpts);

    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
    