% Given a test image, extract s covariance matrices from randomly selected
% regions. Selected regions must have variance above a specific threshold
% to guarantee that they do not contain mainly the background.
% read the image and calculate its covariance descriptor from
% complex wavelet coefficient magnitudes
function [f,xnew] = covcomplex_features7_furkan_v2(I,option)

% extract feature vector of complex wavelet coefficients from a given image
[r c] = size(I);

I = double(I);

if(rem(r,4) ~= 0)
    I = I(1:end-rem(r,4),:);
    r = r - rem(r,4);
end

if(rem(c,4) ~= 0)
    I = I(:,1:end-rem(c,4));
    c = c - rem(c,4);
end

% linear scaling to unit range
x = imnormalize(I,0,1);

% compute the 2-D dual-tree complex wavelet transform of the image
% extract features from 2nd decomposition level because the wavelet is not
% analytic at the first level
Jd = 2;
[Faf, Fsf] = FSfarras;
[af, sf] = dualfilt1;
w = cplxdual2D(x, Jd, Faf, af);

% generate a low-low image the quarter size of the original in each
% dimension to be used in feature vector

% instead you can use 

xnew = imresize(x,1/4);

% simpf = [1/4 1/2 1/4];
% for i = 1:r
%     x(i,:) = conv(x(i,:),simpf,'same');
% end
% 
% for j = 1:c
%     x(:,j) = (conv(x(:,j)',simpf,'same'))';
% end
% 
% xnew = x(2:2:end,2:2:end);
% 
% % two times downsampling and filtering to avoid aliasing
% 
% for i = 1:r/2
%     xnew(i,:) = conv(xnew(i,:),simpf,'same');
% end
% 
% for j = 1:c/2
%     xnew(:,j) = (conv(xnew(:,j)',simpf,'same'))';
% end
% 
% xnew = xnew(2:2:end,2:2:end);
    
if(option == 9)
    
    % intensity, CWT, needle and derivative features
    % 19-dim feature vector
    
    % obtain circular needle output images
    I2 = 255*xnew;
    resI = circularNeedleAnglesSameSize7(I2,10,10);
    
    f = zeros(19,r/4,c/4);
        
    f(1,:,:) = xnew/2;
    
    f(2:9,:,:)=shiftdim(resI(:,:,1:8),2);

    % 2nd level DT-CWT coefficients
    f(10,1:r/4,1:c/4) = sqrt((w{2}{1}{1}{1}(1:r/4,1:c/4)).^2 + (w{2}{2}{1}{1}(1:r/4,1:c/4)).^2);
    f(11,1:r/4,1:c/4) = sqrt((w{2}{1}{1}{2}(1:r/4,1:c/4)).^2 + (w{2}{2}{1}{2}(1:r/4,1:c/4)).^2);
    f(12,1:r/4,1:c/4) = sqrt((w{2}{1}{1}{3}(1:r/4,1:c/4)).^2 + (w{2}{2}{1}{3}(1:r/4,1:c/4)).^2);
    f(13,1:r/4,1:c/4) = sqrt((w{2}{1}{2}{1}(1:r/4,1:c/4)).^2 + (w{2}{2}{2}{1}(1:r/4,1:c/4)).^2);
    f(14,1:r/4,1:c/4) = sqrt((w{2}{1}{2}{2}(1:r/4,1:c/4)).^2 + (w{2}{2}{2}{2}(1:r/4,1:c/4)).^2);
    f(15,1:r/4,1:c/4) = sqrt((w{2}{1}{2}{3}(1:r/4,1:c/4)).^2 + (w{2}{2}{2}{3}(1:r/4,1:c/4)).^2);
   
    
    % norm of 1st and 2nd derivatives
    i = 2:(r/4)-1;
    j = 1:c/4;
    f(16,i,j) = abs(xnew(i+1,j) - xnew(i-1,j))/2;
    f(18,i,j) = abs(2*xnew(i,j) - xnew(i+1,j) - xnew(i-1,j))/2;
    
    j = 2:(c/4)-1;
    i = 1:r/4;
    f(17,i,j) = abs(xnew(i,j+1) - xnew(i,j-1))/2;
    f(19,i,j) = abs(2*xnew(i,j) - xnew(i,j+1) - xnew(i,j-1))/2;
  
    
elseif(option == 10)
    
    % intensity, CWT, needle and derivative features
    % 19-dim feature vector
    
    % obtain circular needle output images
    I2 = 255*xnew;
    resI = circularNeedleMedAnglesSameSize(I2,10);
    
    f = zeros(19,r/4,c/4);
    
    f(1,:,:) = xnew/2;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii);
    end
    
    for i = 1:r/4
        for j = 1:c/4
            
            % 2nd level DT-CWT coefficients
            f(10,i,j) = sqrt(w{2}{1}{1}{1}(i,j)^2 + w{2}{2}{1}{1}(i,j)^2);
            f(11,i,j) = sqrt(w{2}{1}{1}{2}(i,j)^2 + w{2}{2}{1}{2}(i,j)^2);
            f(12,i,j) = sqrt(w{2}{1}{1}{3}(i,j)^2 + w{2}{2}{1}{3}(i,j)^2);
            f(13,i,j) = sqrt(w{2}{1}{2}{1}(i,j)^2 + w{2}{2}{2}{1}(i,j)^2);
            f(14,i,j) = sqrt(w{2}{1}{2}{2}(i,j)^2 + w{2}{2}{2}{2}(i,j)^2);
            f(15,i,j) = sqrt(w{2}{1}{2}{3}(i,j)^2 + w{2}{2}{2}{3}(i,j)^2);
            
            % norm of 1st and 2nd derivatives
            if(i>1 && i<r/4)
                f(16,i,j) = abs(xnew(i+1,j) - xnew(i-1,j))/2;
            end
            if(j>1 && j<c/4)
                f(17,i,j) = abs(xnew(i,j+1) - xnew(i,j-1))/2;
            end
            if(i>1 && i<r/4)
                f(18,i,j) = abs(2*xnew(i,j) - xnew(i+1,j) - xnew(i-1,j))/2;
            end
            if(j>1 && j<c/4)
                f(19,i,j) = abs(2*xnew(i,j) - xnew(i,j+1) - xnew(i,j-1))/2;
            end
            
        end
    end
    
elseif(option == 11)
    
    % intensity, FWT, needle and derivative features
    % 16-dim feature vector
    
    r = r1;
    c = c1;
    
    % obtain circular needle output images
    [I_ll,~,I_sub2,I_sub3_1,I_sub4,~] = directional_FWT2(imstr, 1);
    resI = circularNeedleMaxAnglesSameSize(I_ll,10);
    
    f = zeros(16,r/3,c/3);
    
    f(1,:,:) = I_ll/255;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii);
    end
    
    f(10,:,:) = I_sub2;
    f(11,:,:) = I_sub3_1;
    f(12,:,:) = I_sub4;
    
    for i = 1:r/3
        for j = 1:c/3
            
            % norm of 1st and 2nd derivatives
            if(i>1 && i<r/3)
                f(13,i,j) = abs(I_ll(i+1,j) - I_ll(i-1,j))/2;
            end
            if(j>1 && j<c/3)
                f(14,i,j) = abs(I_ll(i,j+1) - I_ll(i,j-1))/2;
            end
            if(i>1 && i<r/3)
                f(15,i,j) = abs(2*I_ll(i,j) - I_ll(i+1,j) - I_ll(i-1,j))/2;
            end
            if(j>1 && j<c/3)
                f(16,i,j) = abs(2*I_ll(i,j) - I_ll(i,j+1) - I_ll(i,j-1))/2;
            end
            
        end
    end
    
elseif(option == 12)
    
    % intensity, FWT, needle and derivative features
    % 16-dim feature vector
    
    r = r1;
    c = c1;
    
    % obtain circular needle output images
    [I_ll,~,I_sub2,I_sub3_1,I_sub4,~] = directional_FWT2(imstr, 1);
    resI = circularNeedleAnglesSameSize(I_ll,10,10);
    
    f = zeros(16,r/3,c/3);
    
    f(1,:,:) = I_ll/255;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii);
    end
    
    f(10,:,:) = I_sub2;
    f(11,:,:) = I_sub3_1;
    f(12,:,:) = I_sub4;
    
    for i = 1:r/3
        for j = 1:c/3
            
            % norm of 1st and 2nd derivatives
            if(i>1 && i<r/3)
                f(13,i,j) = abs(I_ll(i+1,j) - I_ll(i-1,j))/2;
            end
            if(j>1 && j<c/3)
                f(14,i,j) = abs(I_ll(i,j+1) - I_ll(i,j-1))/2;
            end
            if(i>1 && i<r/3)
                f(15,i,j) = abs(2*I_ll(i,j) - I_ll(i+1,j) - I_ll(i-1,j))/2;
            end
            if(j>1 && j<c/3)
                f(16,i,j) = abs(2*I_ll(i,j) - I_ll(i,j+1) - I_ll(i,j-1))/2;
            end
            
        end
    end
    
elseif(option == 13)
    
    % intensity, FWT, needle and derivative features
    % 16-dim feature vector
    
    r = r1;
    c = c1;
    
    % obtain circular needle output images
    [I_ll,~,I_sub2,I_sub3_1,I_sub4,~] = directional_FWT2(imstr, 1);
    resI = circularNeedleMedAnglesSameSize(I_ll,10);
    
    f = zeros(16,r/3,c/3);
    
    f(1,:,:) = I_ll/255;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii);
    end
    
    f(10,:,:) = I_sub2;
    f(11,:,:) = I_sub3_1;
    f(12,:,:) = I_sub4;
    
    for i = 1:r/3
        for j = 1:c/3
                     
            sum1 = 0;
            if(i>1 && i<r/3 && j>2 && j<r/3-1)
                
                % horizontal
     
                for x = [-1 0 1]
                    for y = [-2 -1 1 2]
                        sum1 = sum1 + I_ll(i+x,j+y);
                    end
                end
                
            end
            
            sum2 = 0;
            if(j>2 && j<r/3-1)
                
                % horizontal
              
                for x = [0]
                    for y = [-2 -1 1 2]
                        sum2 = sum2 + I_ll(i+x,j+y);
                    end
                end
                
            end
            
            sum3 = 0;
            if(j>1 && j<r/3 && i>2 && i<r/3-1)
                
                % vertical
                for y = [-1 0 1]
                    for x = [-2 -1 1 2]
                        sum3 = sum3 + I_ll(i+x,j+y);
                    end
                end
            end
            
            sum4 = 0;
            if(i>2 && i<r/3-1)
                
                % vertical
                for y = [0]
                    for x = [-2 -1 1 2]
                        sum4 = sum4 + I_ll(i+x,j+y);
                    end
                end
            end
            
            % norm of 1st and 2nd derivatives
   
            f(13,i,j) = abs(I_ll(i,j) - sum1/12)/2;
            
            
            f(14,i,j) = abs(I_ll(i,j) - sum2/4)/2;
            
            
            f(15,i,j) = abs(I_ll(i,j) - sum3/12)/2;
            
            
            f(16,i,j) = abs(I_ll(i,j) - sum4/4)/2;

            
        end
    end
    
%     for i = 1:r/3
%         for j = 1:c/3
%             
%             % norm of 1st and 2nd derivatives
%             if(i>1 && i<r/3)
%                 f(13,i,j) = abs(I_ll(i+1,j) - I_ll(i-1,j))/2;
%             end
%             if(j>1 && j<c/3)
%                 f(14,i,j) = abs(I_ll(i,j+1) - I_ll(i,j-1))/2;
%             end
%             if(i>1 && i<r/3)
%                 f(15,i,j) = abs(2*I_ll(i,j) - I_ll(i+1,j) - I_ll(i-1,j))/2;
%             end
%             if(j>1 && j<c/3)
%                 f(16,i,j) = abs(2*I_ll(i,j) - I_ll(i,j+1) - I_ll(i,j-1))/2;
%             end
%             
%         end
%     end


    
elseif(option == 14)
    
    % intensity, 25:1 directional WT with 8 directions, needle
    % 17-dim feature vector
    r = r1;
    c = c1;
    
    % decimate by 5
    I = I2;
    I = I(1:r-rem(r,5),1:c-rem(c,5));
    r = r - rem(r,5);
    c = c - rem(c,5);
    
    d = 5;
    Ir = zeros(r,c/d);
    
    for i = 1:r
        Ir(i,:) = decimate(I(i,:),d);
    end
    
    I_decimated = zeros(r/d,c/d);
    for j = 1:c/d
        I_decimated(:,j) = decimate(Ir(:,j),d);
    end
    
    I_decimated = I_decimated/255;
    I = I/255;
    
    % obtain 25:1 subbands and circular needle output images
    subbands = directionalWT5(I);
    resI = circularNeedleMedAnglesSameSize(I_decimated,10);
    
    f = zeros(17,r/5,c/5);
    
    f(1,:,:) = I_decimated*0.784;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii);
        f(ii+9,:,:) = subbands{ii};
    end
    
    xnew = f(1,:,:);
    
elseif(option == 15)
    
    % intensity, 25:1 directional WT with 8 directions, needle
    % 17-dim feature vector
    r = r1;
    c = c1;
    
    % decimate by 5
    I = I2;
    I = I(1:r-rem(r,5),1:c-rem(c,5));
    r = r - rem(r,5);
    c = c - rem(c,5);
    
    d = 5;
    Ir = zeros(r,c/d);
    
    for i = 1:r
        Ir(i,:) = decimate(I(i,:),d);
    end
    
    I_decimated = zeros(r/d,c/d);
    for j = 1:c/d
        I_decimated(:,j) = decimate(Ir(:,j),d);
    end
    
    I_decimated = I_decimated/255;
    I = I/255;
    
    % obtain 25:1 subbands and circular needle output images
    subbands = directionalWT5(I);
    resI = circularNeedleMaxAnglesSameSize(I_decimated,10);
    
    f = zeros(17,r/5,c/5);
    
    f(1,:,:) = I_decimated*0.784;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii);
        f(ii+9,:,:) = subbands{ii};
    end
    
    xnew = f(1,:,:);
    
elseif(option == 16)
    
    % intensity, 25:1 directional WT with 8 directions, needle
    % 17-dim feature vector
    r = r1;
    c = c1;
    
    % decimate by 5
    I = I2;
    I = I(1:r-rem(r,5),1:c-rem(c,5));
    r = r - rem(r,5);
    c = c - rem(c,5);
    
    d = 5;
    Ir = zeros(r,c/d);
    
    for i = 1:r
        Ir(i,:) = decimate(I(i,:),d);
    end
    
    I_decimated = zeros(r/d,c/d);
    for j = 1:c/d
        I_decimated(:,j) = decimate(Ir(:,j),d);
    end
    
    
    I = I/255;
    
    % obtain 25:1 subbands and circular needle output images
    subbands = directionalWT5(I);
    resI = circularNeedleAnglesSameSize(I_decimated,10,10);
    
    I_decimated = I_decimated/255;
    
    f = zeros(17,r/5,c/5);
    
    f(1,:,:) = I_decimated*0.784;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii)/50;
        f(ii+9,:,:) = subbands{ii};
    end
    
    xnew = f(1,:,:);
    
elseif(option == 17)
    
    % intensity, our ST-CWT, needle and derivative features
    % 19-dim feature vector
    
    xnew = decimate_micro(x2,16);
    [rx cx] = size(xnew);
    
    % obtain circular needle output images
    I2 = 255*xnew;
    resI = circularNeedleAnglesSameSize(I2,10,10);
    
    % feature tensor
    f = zeros(19,rx,cx);
    
%     size(resI)
%     size(x2)
%     size(xnew)
    
    % use the 16-by-16 decimated image
    f(1,:,:) = xnew/2;
    
    % needle output image pixel values
    for ii = 1:8
        f(ii+1,:,:) = resI(:,:,ii);
    end
    
    % perform two-level 2-D ST-CWT decomposition
    subbands = balanced_lifting_switched_2d_2ndlevel(x2);
    
    for i = 1:rx
        for j = 1:cx
            
            % 2nd level ST-CWT coefficients
            f(10,i,j) = sqrt(subbands{2,1}(i,j)^2 + subbands{3,3}(i,j)^2);
            f(11,i,j) = sqrt(subbands{2,2}(i,j)^2 + subbands{3,4}(i,j)^2);
            f(12,i,j) = sqrt(subbands{2,3}(i,j)^2 + subbands{4,1}(i,j)^2);
            f(13,i,j) = sqrt(subbands{2,4}(i,j)^2 + subbands{4,2}(i,j)^2);
            f(14,i,j) = sqrt(subbands{3,1}(i,j)^2 + subbands{4,3}(i,j)^2);
            f(15,i,j) = sqrt(subbands{4,4}(i,j)^2 + subbands{3,2}(i,j)^2);
            
            % norm of 1st and 2nd derivatives
            if(i>1 && i<rx)
                f(16,i,j) = abs(xnew(i+1,j) - xnew(i-1,j))/2;
            end
            if(j>1 && j<cx)
                f(17,i,j) = abs(xnew(i,j+1) - xnew(i,j-1))/2;
            end
            if(i>1 && i<rx)
                f(18,i,j) = abs(2*xnew(i,j) - xnew(i+1,j) - xnew(i-1,j))/2;
            end
            if(j>1 && j<cx)
                f(19,i,j) = abs(2*xnew(i,j) - xnew(i,j+1) - xnew(i,j-1))/2;
            end
            
        end
    end
    
end

% % construct a feature vector for each pixel of the image using the
% % magnitudes of the complex wavelet coefficients in six directional
% % subbands and the pixel intensity and x and y derivatives
% f = zeros(5,r/2,c/2);
%
% % the approximation image is obtained by averaging the four approximation
% % subbands obtained by dual-tree CWT
% Ilow = (1/4)*(w{2}{2}{2}+w{2}{2}{1}+w{2}{1}{2}+w{2}{1}{1});





