function [H2, H4, H1] = calculateFilterBanks(img, type, numberOfStatistics)
% This code is copyrighted by Stanford University.
% Author: Ashutosh Saxena, Justin Drieymeyer
% Code cannot be used for commercial purposes.
% Use of the code, if used with permission from authors, must acknowledge:
% Learning to Grasp Novel Objects using Vision,
% Ashutosh Saxena, Justin Driemeyer, Justin Kearns, Chioma Osondu, Andrew Y. Ng. In 10th International Symposium on Experimental Robotics, ISER, 2006.
% Robotic Grasping of Novel Obects,
% Ashutosh Saxena, Justing Driemeyer, Justin Kearns, Andrew Y. Ng, In Neural Information Processing Systems (NIPS), 2006.
% More details at: http://ai.stanford.edu/~asaxena/learninggrasp/
% Version 0.1:   Aug 2006.


if nargin < 3
    numberOfStatistics = 1;
end

persistent L3 E3 S3 gMask		% the Laws' masks
persistent NB1 NB2 NB3 NB4 NB5 NB6

global nLaw 

%nLaw                 = 0;
%nRadon               = 0;
%nHarris              = 0;
%harrisSmoothingWidth = 9;

switch lower(type)
    case 'purelaws'
        nLaw = 9;
    case {'laws', 'texture'}
        nLaw = 17;
    case {'harris', 'texture'}
        nHarris = 10;
    case {'radon'}
        nRadon = nTopRadon*180/thetaStep;
    case {'lawsradon'}
        nLaw = 17;
        nRadon = nTopRadon*180/thetaStep;
    case {'lawsharris'}
        nLaw = 17;
        nHarris = 10;
end
%nDim = nLaw + nHarris;
 
%convultional filters
if length(L3) <= 0
    MatrixConditioner = .2;
    L3 = [1 2 1]/128/MatrixConditioner;      
    E3 = [-1 0 1];
    S3 = [-1 2 -1];
    %L5 = [1 2 3 2 1];       %can be improved
    %L7 = [1 2 3 4 3 2 1];   % can be improved
    gMask = gaussMask(4,1.4);
    
    %Navatia Babu filters oriented at six angles (0,30,60,90,120,150)
    NB1 = [ -100, -100, 0, 100, 100; ...
            -100, -100, 0, 100, 100; ...
            -100, -100, 0, 100, 100; ...
            -100, -100, 0, 100, 100; ...
            -100, -100, 0, 100, 100];
    NB2 = [ -100, 32, 100, 100, 100; ...
            -100,-78, 92,  100, 100; ...
            -100,-100, 0,  100, 100; ...
            -100,-100,-92, 78,  100; ...
            -100,-100,-100,-32, 100];
    
    NB1 = NB1/2000;
    NB2 = NB2/2000;
        
    NB3 = -NB2';
    NB4 = -NB1';
    NB5 = -NB3(end:-1:1,:);
    NB6 = NB5'; 
end

ycbcr_img = double(rgb2ycbcr(img)); %This computer has no image processing functions!!

%ycbcr_img(:,:,1) = ycbcr_img(:,:,1) - mean( mean( ycbcr_img(:,:,1) ));
%ycbcr_img(:,:,2) = ycbcr_img(:,:,2) - mean( mean( ycbcr_img(:,:,2) ));
%ycbcr_img(:,:,3) = ycbcr_img(:,:,3) - mean( mean( ycbcr_img(:,:,3) ));

imgY = ycbcr_img(:,:,1);% - mean( mean( ycbcr_img(:,:,1) ));
%imgEdge = edge(imgY, 'sobel');

clear img;
H = zeros( size(imgY,1), size(imgY,2), nLaw);

%L3img = conv2(imgY, L3, 'same');
%E3img = conv2(imgY, E3, 'same');
%S3img = conv2(imgY, S3, 'same');
%E3imgT = conv2(imgY, E3', 'same');
   
H(:,:,1) = conv2(imgY, L3'*L3, 'same');
H(:,:,2) = conv2(imgY, L3'*E3, 'same');
H(:,:,3) = conv2(imgY, L3'*S3, 'same');
   
H(:,:,4) = conv2(imgY, E3'*L3, 'same');
H(:,:,5) = conv2(imgY, E3'*E3, 'same');
H(:,:,6) = conv2(imgY, E3'*S3, 'same');
   
H(:,:,7) = conv2(imgY, S3'*L3, 'same');
H(:,:,8) = conv2(imgY, S3'*E3, 'same');
H(:,:,9) = conv2(imgY, S3'*S3, 'same');
   
H(:,:,10) = conv2( (ycbcr_img(:,:,2)), L3'*L3, 'same');
H(:,:,11) = conv2( (ycbcr_img(:,:,3)), L3'*L3, 'same');

%H(:,:,12) = conv2(E3img, L5', 'same');
%H(:,:,13) = conv2(E3img, L7', 'same');
%H(:,:,13) = conv2(E3imgT, L5, 'same');
%H(:,:,15) = conv2(E3imgT, L7, 'same');

H(:,:,12) = conv2(imgY, NB1, 'same');
H(:,:,13) = conv2(imgY, NB2, 'same');
H(:,:,14) = conv2(imgY, NB3, 'same');
H(:,:,15) = conv2(imgY, NB4, 'same');
H(:,:,16) = conv2(imgY, NB5, 'same');
H(:,:,17) = conv2(imgY, NB6, 'same');


H1 = 0;%abs(H);
H2 = H.^2;      %Energy
H4 = 0;%H2.^2;      %Kurtosis

%=======for Harris texture gradient ===========
%Ix2 = conv2(H2(:,:,12), gMask, 'same');
%Iy2 = conv2(H2(:,:,14), gMask, 'same');
%Ixy = conv2(H(:,:,12) .* H(:,:,14), gMask, 'same');

%anglesHarris = zeros( size(Ix2,1), size(Ix2,2), 2);
%eigValue = zeros( size(Ix2,1), size(Ix2,2), 2);
%upsilon = 1e-23;

%gradientMagImage = sqrt(H2(:,:,12) + H2(:,:,14) );
%angleImage = round( (nHarris-1) * (pi/2+atan( H(:,:,12) ./ (H(:,:,14)+upsilon))) /pi) + 1;