%Description: G Filters
%Usage:
%       [Gx,Gy]=gfilter(lambda,tau,filterLength)
%       lambda>0
%       0.0 <= tau <= 1.0
%       filterLength >= 1  
function [Gx,Gy]=gfilter(lambda,tau,filterLength)

if tau > 1.0 
    disp('Tau must be less than or equal to one') ;
    tau=1.0 ;
end
if tau < 0.0
    disp('Tau must be greater than or equal to zero') ;    
    tau=0.0 ;
end    
if lambda <= 0.0
    disp('Lambda must be greater than zero') ;    
    return ;
end        
if filterLength <= 1 
    disp('Filter Length must be greater than one') ;
    return ;
end

a     = lambda * tau ;
b     = lambda * (1.0 - tau) ;
delta = b*b - 4.0 * a ; 
% fprintf('\na=%f,b=%f,delta=%f',a,b,delta) ;

Gx     = zeros(2*filterLength+1,2*filterLength+1) ;
Gy     = zeros(2*filterLength+1,2*filterLength+1) ;

if a == 0 
%     fprintf('\nCase: A==0') ;                            
    alpha = 1.0 / sqrt(2.0*b) ;
%     fprintf('\nalpha=%f',alpha) ;                    
    for m = -filterLength:filterLength
        for n = -filterLength:filterLength        
            k = abs(m) +  abs(n) ;
            Gy(m+filterLength+1,n+filterLength+1) = -sign(m)*exp(-alpha*k) ;
            Gx(m+filterLength+1,n+filterLength+1) = -sign(n)*exp(-alpha*k) ;
        end
    end                
%     sumOfGx= sum(sum(Gx)) ;
%     sumOfGy= sum(sum(Gy)) ;
%     fprintf('\nSum of the Gx filter: %f',sumOfGx);
%     fprintf('\nSum of the Gy filter: %f',sumOfGy);    
elseif b == 0 
%     fprintf('\nCase: b=0') ;                
    c = 1.0 / sqrt(sqrt(4.0*a)) ;
    theta = pi/4.0 ;
    phi = sin(theta) *c ;
    alpha = cos(theta)*c ;

%     fprintf('\ntheta=%f \nphi=%f \nalpha=%f',theta,phi,alpha) ;        
    for m = -filterLength:filterLength
        for n = -filterLength:filterLength        
            k = abs(m) +  abs(n) ;
            Gy(m+filterLength+1,n+filterLength+1) = -c * sign(m) * ((cos(2*theta)*cos(phi*k))+(sin(2*theta)*sin(phi*k)))*exp(-alpha*k) ;
            Gx(m+filterLength+1,n+filterLength+1) = -c * sign(n) * ((cos(2*theta)*cos(phi*k))+(sin(2*theta)*sin(phi*k)))*exp(-alpha*k) ;
        end
    end              
%     sumOfGx= sum(sum(Gx)) ;
%     sumOfGy= sum(sum(Gy)) ;
%     fprintf('\nSum of the Gx filter: %f',sumOfGx);
%     fprintf('\nSum of the Gy filter: %f',sumOfGy);            
else
if delta > 0
%     fprintf('\nCase: Delta>0') ;        
    e1 = sqrt((b+sqrt(delta))/(4.0*a)) ;
    e2 = sqrt((b-sqrt(delta))/(4.0*a)) ;
%     fprintf('\ne1=%f,e2=%f',e1,e2) ;    
    for m = -filterLength:filterLength
        for n = -filterLength:filterLength        
            k = abs(m) +  abs(n) ;
            Gy(m+filterLength+1,n+filterLength+1) = sign(m) * (exp(-e1*k) - exp(-e2*k)) ;
            Gx(m+filterLength+1,n+filterLength+1) = sign(n) * (exp(-e1*k) - exp(-e2*k)) ;
        end
    end    
%     sumOfGx= sum(sum(Gx)) ;
%     sumOfGy= sum(sum(Gy)) ;
%     fprintf('\nSum of the Gx filter: %f',sumOfGx);
%     fprintf('\nSum of the Gy filter: %f',sumOfGy);            
else
if delta < 0
%     fprintf('\nCase: Delta<0') ;                
    c = 1.0 / sqrt(sqrt(4.0*a)) ;
    theta = 0.5 * atan(sqrt(abs(delta))/b) ;
    phi = sin(theta) *c ;
    alpha = cos(theta)*c ;

%     fprintf('\ntheta=%f \nphi=%f \nalpha=%f',theta,phi,alpha) ;        
    for m = -filterLength:filterLength
        for n = -filterLength:filterLength        
            k = abs(m) +  abs(n) ;
            Gy(m+filterLength+1,n+filterLength+1) = -c * sign(m) * ((cos(2*theta)*cos(phi*k))+(sin(2*theta)*sin(phi*k)))*exp(-alpha*k) ;
            Gx(m+filterLength+1,n+filterLength+1) = -c * sign(n) * ((cos(2*theta)*cos(phi*k))+(sin(2*theta)*sin(phi*k)))*exp(-alpha*k) ;
        end
    end              
%     sumOfGx= sum(sum(Gx)) ;
%     sumOfGy= sum(sum(Gy)) ;
%     fprintf('\nSum of the Gx filter: %f',sumOfGx);
%     fprintf('\nSum of the Gy filter: %f',sumOfGy);            
else
if delta == 0
%     fprintf('\nCase: Delta==0') ;                            
    alpha = 1.0/ sqrt(b) ;
    for m = -filterLength:filterLength
        for n = -filterLength:filterLength        
            k = abs(m) +  abs(n) ;
            Gy(m+filterLength+1,n+filterLength+1) = -sign(m)*(alpha*abs(n)+1.0) * exp(-alpha*k) ;
            Gx(m+filterLength+1,n+filterLength+1) = -sign(n)*(alpha*abs(m)+1.0) * exp(-alpha*k) ;
        end
    end    
%     sumOfGx= sum(sum(Gx)) ;
%     sumOfGy= sum(sum(Gy)) ;
%     fprintf('\nSum of the Gx filter: %f',sumOfGx);
%     fprintf('\nSum of the Gy filter: %f',sumOfGy);                            
            end
        end
    end
end