% h = ss_DoGauss( k, scale, sigma, theta )
% ------------------------------------
%
% Component of Scale Space Edges, version 1.0.
%
% Generates a single Difference of Gauss (DoG) filter.
%
% Input:
%
%     k		- [0, 3(4)] kth order derivative. 
%
%     scale	- [2,4] Separation.
%
%     sigma	- Standard deviation.
%
%     theta	- Rotation angle of filter, ccw (rad). 
%
% Output:
%
%    Returns a set of triples, containing the 
%    coordinates and values of a discrete sampling of the kth 
%    derivative in one direction of a 2D gaussian with 
%    a standard deviation of sigma*scale.  The array is normalized 
%    for a maximal response of .5 (to 1 for full map).  Half the map 
%    (positive coordinates, and half value at the odd point of 0,0)
%    is generated, as these are all bilaterally symmetric or 
%    anti-symmetric. 
%
%
% Mark Dow, last modified January 2002
% University of Oregon
% Brain Development Lab
% dow@braindev.uoregon.edu
%

function h = ss_DoGauss( k, scale, sigma, theta );



%max_width = 5;
max_width = 6;
% NOTE: If k = 4, max_width should be at least 6


%n = 2*round( (2*sigma*max_width + 1)/2 ) + 1;
n = 25;

if k >= 4

    'WARNING: 4th derivative is not balanced.'
    'The second moment should equal zero ( sum(hi*ri^2) = 0 )?'
end
if ~( k == 0  |  k == 1  |  k == 2  |  k == 3  |  k == 4 )
    a = 'Currently can deal with only 0-4th derivatives.'
    h = 0;
    return;
    
end

r = [ cos(theta) -sin(theta);
      sin(theta)  cos(theta) ];

m = 0;  
 
for i = 0 : (n+1)/2
    for j = - (n+1)/2 : (n+1)/2
        
        if ~( i == 0  &  j < 0 ) ...            % Half of center line.
        
            fac = 1;                    
            if i == 0 & j == 0                  % Center pixel has half weight.
                fac = .5; 
            end
        
            m = m + 1;       
        
            u = r * [ j i ]';
        
            h( m, 1 ) = i;
            h( m, 2 ) = j;
        
            if k == 0
                h( m, 3 ) = fac*gauss( u(1), sigma )*gauss( u(2), sigma );
            end 
            if k == 1
%                h( m, 3 ) = fac*gauss( u(1), scale*sigma )*dgauss( u(2), scale*sigma );
                h( m, 3 ) = fac*gauss( u(1), sigma )*( gauss( u(2) - scale/2, sigma ) ...
                                                     - gauss( u(2) + scale/2, sigma ) );
            end    
            if k == 2
%                h( m, 3 ) = fac*gauss( u(1), scale*sigma )*d2gauss( u(2), scale*sigma );
                h( m, 3 ) = fac*gauss( u(1), sigma )*( -   gauss( u(2) - scale, sigma ) ...
                                                       + 2*gauss( u(2), sigma ) ...
                                                       -   gauss( u(2) + scale, sigma ) );
            end
            if k == 3
%                h( m, 3 ) = fac*gauss( u(1), scale*sigma )*d3gauss( u(2), scale*sigma );
                h( m, 3 ) = fac*gauss( u(1), sigma )*( +   gauss( u(2) - 3*scale/2, sigma ) ...
                                                       - 3*gauss( u(2) - scale/2, sigma ) ...
                                                       + 3*gauss( u(2) + scale/2, sigma ) ...
                                                       -   gauss( u(2) + 3*scale/2, sigma ));
            end
            if k == 4
                h( m, 3 ) = fac*gauss( u(1), scale*sigma )*d4gauss( u(2), scale*sigma );
            end
        end
    end
end


% Balancing:
if k == 3

    % Net moment = 0.
    moment = sum( h(:,3).*( h(:,1)*cos(theta) ) );
    h(:,3) = h(:,3) - ( sign(h(:,1)).*abs( h(:,3) )*moment )/( sum( sign(h(:,1)).*abs( h(:,3) ).*( h(:,1)*cos(theta) ) ) );
end

% Normalization:
if k == 0  |  k == 1  |  k == 3

    % Normalize to max response of .5.
    h(:,3) = h(:,3) ./ (2*sum(abs(h(:,3))));
end
if k == 2

    % Normalize positive values to .25 and negative values to .25.
    i_pos = find( h(:,3) > 0 );
    i_neg = find( h(:,3) < 0 );
    h(i_pos,3) = h(i_pos,3) ./ ( 4*sum( h(i_pos,3) ));
    h(i_neg,3) = h(i_neg,3) ./ (-4*sum( h(i_neg,3) ));
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function "gauss.m":

function y = gauss(x,std)

%y = exp(-x^2/(2*std^2)) / (std*sqrt(2*pi));
y = exp(-x^2/(2*std^2));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function "dgauss.m" (first order derivative):

function y = dgauss(x,std)

%y = x * gauss(x,std) / std^2;
y = x * gauss(x,std);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function "d2gauss.m" (second order derivative):

function y = d2gauss(x,std)

y = -(x^2 - std^2) * gauss(x,std);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function "d3gauss.m" (third order derivative):

function y = d3gauss(x,std)

y = -(3*std^2*x - x^3 ) * gauss(x,std);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function "d4gauss.m" (fourth order derivative):

function y = d4gauss(x,std)

y = ( x^4 - 6*x^2*std^2 + 3*std^4 ) * gauss(x,std);