% 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 newim = ridgefilter(im, orient, freq, kx, ky, showfilter)

    if nargin == 5
        showfilter = 0;
    end
    
    angleInc = 3;  
    
    im = double(im);
    [rows, cols] = size(im);
    newim = zeros(rows,cols);
    
    [validr,validc] = find(freq > 0);  
    ind = sub2ind([rows,cols], validr, validc);

    freq(ind) = round(freq(ind)*100)/100;
    
    unfreq = unique(freq(ind)); 
    
    freqindex = ones(100,1);
    for k = 1:length(unfreq)
        freqindex(round(unfreq(k)*100)) = k;
    end
    
    filter = cell(length(unfreq),180/angleInc);
    sze = zeros(length(unfreq),1);
    
    for k = 1:length(unfreq)
        sigmax = 1/unfreq(k)*kx;
        sigmay = 1/unfreq(k)*ky;
        
        sze(k) = round(3*max(sigmax,sigmay));
        [x,y] = meshgrid(-sze(k):sze(k));
        reffilter = exp(-(x.^2/sigmax^2 + y.^2/sigmay^2)/2)...
                .*cos(2*pi*unfreq(k)*x);

        for o = 1:180/angleInc
            filter{k,o} = imrotate(reffilter,-(o*angleInc+90),'bilinear','crop'); 
        end
    end

    if showfilter 
        figure(7), imshow(filter{1,end},[]); title('filter'); 
    end
    
    maxsze = sze(1);    
    finalind = find(validr>maxsze & validr<rows-maxsze & ...
                    validc>maxsze & validc<cols-maxsze);
    
    maxorientindex = round(180/angleInc);
    orientindex = round(orient/pi*180/angleInc);
    i = find(orientindex < 1);   orientindex(i) = orientindex(i)+maxorientindex;
    i = find(orientindex > maxorientindex); 
    orientindex(i) = orientindex(i)-maxorientindex; 

    for k = 1:length(finalind)
        r = validr(finalind(k));
        c = validc(finalind(k));

        filterindex = freqindex(round(freq(r,c)*100));
        
        s = sze(filterindex);   
        newim(r,c) = sum(sum(im(r-s:r+s, c-s:c+s).*filter{filterindex,orientindex(r,c)}));
    end

    
