% ORIENTATIONFILTER  Generate orientation selective filterings of an image
%
% Usage:  oim = orientationfilter(im, norient, angoverlap, boost, cutoff, histcut)
%
% Arguments:     im - Image to be filtered.
%           norient - Number of orientations, try 8.
%        angoverlap - Angular bandwidth overlap factor.  A value of 1 gives a 
%                     minimum overlap between adjacent filter orientations
%                     while still providing full coverage of the spectrum.
%                     However, it is often useful to provide greater overlap
%                     between adjacent filter orientations, especially if the
%                     output is to be viewed using the CYCLEMIX image
%                     blender.  This gives smoother transitions between
%                     images.  Try values 1.5 to 2.0 
%             boost - The ratio that high frequency values are boosted
%                     (or supressed) relative to the low frequency values.
%                     Try values in the range 2 - 4 to start with. Use a
%                     value of 1 for no boost. 
%            cutoff - Cutoff frequency of the highboost filter 0 - 0.5, try 0.2
%           histcut - Percentage of the histogram extremes to truncate.  This
%                     is useful in preventing outlying values in the data
%                     from dominating in the image rescaling for display. Try
%                     a small amount, say, 0.01%
%
% Returns:      oim - Cell array of length 'norient' containing the
%                     orientation filtered images.
%
%
% This function is designed to help identify oriented structures within an image
% by applying orientation selective filters.  If one is looking for lineaments
% it is typically useful to boost the high frequency components in the image as
% well, hence the combination of the two filters
%
% Example: Generate 8 orientation filterings of an image with an angular
% bandwidth overlap factor of 1.5 and amplifying spatial frequencies greater
% than 0.1 (wavelenths < 10 pixels) by a factor of 2.  Also truncate the image
% histograms so that 0.01% of image pixels are saturated.  This helps ensure
% that outlying data values do not dominate when the image is scaled for
% display.
%
% >> oim = orienationfilter(im, 8, 1.5, 2, 0.1, 0.01);
%
% An effective way to view all the output images is to use the cyclic image
% blending tool CYCLEMIX.m
%
% >> cyclemix(oim);
%
% Note that the cyclemix control wheel varies from 0 - 2pi and these angles are
% mapped to the orientation filtered images which vary from 0 - pi in angle.
% This makes the interface a bit counterintuitive to use.  To get around this
% one can replicate the the cell array of orientation filtered images giving an
% array that corresponds to the angles 0 - pi, 0 - pi.  When this is passed to
% cyclemix you get a much more intuitive interface
%
% >> cyclemix({oim{:} oim{:}})
% 
% See also: CYCLEMIX, HIGHBOOSTFILTER, HISTTRUNCATE 

% Copyright (c) 2012-2017 Peter Kovesi
% www.peterkovesi.com
% 
% Permission is hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to deal
% in the Software without restriction, subject to the following conditions:
% 
% The above copyright notice and this permission notice shall be included in 
% all copies or substantial portions of the Software.
%
% The Software is provided "as is", without warranty of any kind.

% August 2012   Original version
% July   2014   Revised to incorporate highboost filtering and histogram
%               truncation.
% December 2017 Eliminated use of fliplr etc to create symmetric filter as
%               this was slightly incorrect.  Also ensured that filters have
%               0 DC value. 


function oim = orientationfilter(im, norient, angoverlap, boost, cutoff, histcut)
    
    IM = fft2(im);
    [rows,cols] = size(im);

    if ~exist('angoverlap','var')
        angScale = norient/2;
    else
        angScale = norient/2 / angoverlap;
    end
    
    % If a highboost filter has been specified apply it to the fourier
    % transform of the image.
    if exist('cutoff', 'var')
        IM = IM .* highboostfilter([rows,cols], cutoff, 1, boost);
    end
    
    if ~exist('histcut','var'), histcut = 0; end
    
    % Construct frequency domain filter matrices
    [ ~ , u1, u2] = filtergrid(rows,cols);
    
    theta = atan2(-u2,u1);            % Matrix values contain polar angle.
                                      % (note -ve y is used to give +ve
                                      % anti-clockwise angles)
    sintheta = sin(theta);
    costheta = cos(theta);

    for o = 1:norient                 % For each orientation...
        angl = (o-1)*pi/norient;           % Filter angle.
        
        % For each point in the filter matrix calculate the angular distance from
        % the specified filter orientation.  To overcome the angular wrap-around
        % problem sine difference and cosine difference values are first computed
        % and then the atan2 function is used to determine angular distance.
        ds = sintheta * cos(angl) - costheta * sin(angl);    % Difference in sine.
        dc = costheta * cos(angl) + sintheta * sin(angl);    % Difference in cosine.
        dtheta = abs(atan2(ds,dc));                          % Absolute angular distance.
      
        % Make the filter symmetric so that the angle deviation from angl =
        % deviation from angle + pi.
        dtheta(dtheta>pi/2) = pi - dtheta(dtheta>pi/2);
    
        % Scale theta so that cosine spread function has the right wavelength
        % and clamp to pi
        dtheta = min(dtheta*angScale, pi);

        % The orientation filter haa a cosine cross section in the frequency domain
        % The spread function is cos(dtheta) between -pi and pi.  We add 1,
        % and then divide by 2 so that the value ranges 0-1
        spread = (cos(dtheta)+1)/2;        
        spread(1,1) = 0;  % Ensure 0 DC
        % show(fftshift(spread),o);
        
        % Apply orientation filter
        oim{o} = real(ifft2(IM.*spread));
        
        % Truncate extremes of image histogram
        if histcut
            oim{o} = histtruncate(oim{o}, histcut, histcut);
        end
    end