function [featureVector] = makeFeatureTexture(img, type, scale)
%this function returnd and displays feature vector from the image 
%startMakingNewFeatureTexture should be executed once before calling this
%function

% Version 0.1, Feb 2006.
% This code is copyrighted by Ashutosh Saxena, Stanford University
% It is available for non-commercial use only. For commericial use, please
% contact the authors.
% Any publication or report resulting from Use of this code, should cite:
% Learning Depth from Single Monocular Images, Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng, NIPS 2005.
% 3-D Depth Reconstruction from a Single Still Image, Ashutosh Saxena, Sung H. Chung, Andrew Y. Ng, To appear in IJCV 2007.
% For more info and dataset, visit: http://ai.stanford.edu/~asaxena/learningdepth/
% For a better version of the algorithm, also visit: 
% http://ai.stanford.edu/~asaxena/reconstruction3d/



persistent L3 E3 S3 L5 L7 gMask       % the Laws' masks

%global nTopRadon thetaStep
global columnWidth nGridCol nGridRow rowWidth
global scale2Factor

if nargin < 3
    scale = 1;
end
if nargin  < 2
    type = 'texture';
end

if scale == 1
    cWidth = columnWidth;
    rWidth = rowWidth;
    nGCol = nGridCol;
    nGRow = nGridRow;
elseif scale == 2
    cWidth = columnWidth / scale2Factor;
    rWidth = rowWidth / scale2Factor;
    nGCol = nGridCol / scale2Factor;
    nGRow = nGridRow / scale2Factor;
end

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

switch lower(type)
    case 'purelaws'
        nLaw = 9;
    case {'laws', 'texture'}
        nLaw = 15;
    case {'harris', 'texture'}
        nHarris = 10;
    case {'radon'}
        nRadon = nTopRadon*180/thetaStep;
    case {'lawsradon'}         % 'texture' means the best one, which right now is lawsradon
        nLaw = 15;
        nRadon = nTopRadon*180/thetaStep;
    case {'lawsharris'}
        nLaw = 15;
        nHarris = 10;
end
    
%convultional filters
if length(L3) <= 0
    L3 = [1 2 1];
    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);
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');

H = zeros( size(imgY,1), size(imgY,2), 15);

L3img = conv2(imgY, L3, 'same');
E3img = conv2(imgY, E3, 'same');
S3img = conv2(imgY, S3, 'same');
E3imgT = conv2(imgY, E3', 'same');
   
H(:,:,1) = conv2(L3img, L3', 'same');
H(:,:,2) = conv2(E3img, L3', 'same');
H(:,:,3) = conv2(S3img, L3', 'same');
   
H(:,:,4) = conv2(L3img, E3', 'same');
H(:,:,5) = conv2(E3img, E3', 'same');
H(:,:,6) = conv2(S3img, E3', 'same');
   
H(:,:,7) = conv2(L3img, S3', 'same');
H(:,:,8) = conv2(E3img, S3', 'same');
H(:,:,9) = conv2(S3img, 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(:,:,14) = conv2(E3imgT, L5, 'same');
H(:,:,15) = conv2(E3imgT, L7, 'same');
  
H2 = H.^2;      %Energy
H4 = 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;

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


%forming spatial masks



%=========================================
%featureVectorTmp = zeros( nGRow, nGCol, 2*nLaw+3*nHarris);
featureVectorTmp = zeros( nGRow, nGCol, 2*nLaw+nHarris);

for g = 1:nGRow
    for c = 1:nGCol
        gridLeft = round( g*rWidth-(rWidth-1) );
        gridRight = round(g*rWidth);
        gridTop = round(c*cWidth-(cWidth-1));
        gridBot = round(c*cWidth);
        
        %LAWS-----------
        if nLaw
            featureVectorTmp(g,c,1:nLaw) = sum( sum( H2( gridLeft:gridRight, gridTop:gridBot, 1:nLaw), 2), 1);
            featureVectorTmp(g,c,nLaw+(1:nLaw)) = sum( sum( H4( gridLeft:gridRight, gridTop:gridBot, 1:nLaw), 2), 1);
        end
        if nHarris
            %Can be optimized
            for row = gridLeft:gridRight
                for col = gridTop:gridBot;
                    featureVectorTmp(g,c,2*nLaw+angleImage(row,col) ) = featureVectorTmp(g,c,2*nLaw+angleImage(row,col) ) + ...
                                            gradientMagImage(row,col);
                    %featureVectorTmp(g,c,2*nLaw+nHarris+angleImage(row,col) ) = featureVectorTmp(g,c,2*nLaw+nHarris+angleImage(row,col) ) + ...
                    %                        gradientMagImage(row,col) .* (2*col<(gridTop+gridBot));
                    %featureVectorTmp(g,c,2*nLaw+2*nHarris+angleImage(row,col) ) = featureVectorTmp(g,c,2*nLaw+2*nHarris+angleImage(row,col) ) + ...
                    %                        gradientMagImage(row,col) .* (2*col>=(gridTop+gridBot));
                                        
                end
            end
            featureVectorTmp(g,c,2*nLaw+(1:nHarris) ) = featureVectorTmp(g,c,2*nLaw+(1:nHarris) ) ./ ...
                (max( featureVectorTmp(g,c,2*nLaw+(1:nHarris) ) ) + upsilon);
            %featureVectorTmp(g,c,2*nLaw+nHarris+(1:nHarris) ) = featureVectorTmp(g,c,2*nLaw+nHarris+(1:nHarris) ) ./ ...
            %    (max( featureVectorTmp(g,c,2*nLaw+nHarris+(1:nHarris) ) ) + upsilon);
            %featureVectorTmp(g,c,2*nLaw+2*nHarris+(1:nHarris) ) = featureVectorTmp(g,c,2*nLaw+2*nHarris+(1:nHarris) ) ./ ...
            %    (max( featureVectorTmp(g,c,2*nLaw+2*nHarris+(1:nHarris) ) ) + upsilon);
        end
    end
end
   
featureVector = featureVectorTmp;   %permute( featureVectorTmp, [3 1 2]);
%featureVector = reshape( featureVector3, [size(featureVector2,1)*size(featureVector3,2) nGridCol] );

%    if displayFlag
%        figure, 
%        subplot(2,2,1), imshow(img);
%        subplot(2,2,2), imagesc(ycbcr_img(:,:,1));
%        subplot(2,2,4), imagesc(ycbcr_img(:,:,2));
%        subplot(2,2,3); plot( squeeze(featureVector(frame,:,:)) );
%        axis tight; 
%    end
   
return;