% FilterConv4.m
% Version of 30 June 2021
% Gabor Filter Load and Convolve with Image for VisNet. ETR Oct 2014 from FilterConv3.m of April 2012.
% e.g. FilterConv4('ramp.dat','',0); % standard for use with VisNet 2021
% ramp.dat is created with CreateRawImage.m

% e.g. FilterConv4('image256.dat','/home/erolls/imagedata/tmp/',0)
% Loads the complete set of gabor filters currently output by MakeGaborFilters.m
% (assumes the filter files are in the current path)
% For VisNet, ang=0 is horizontal
% Visnet images are uint8 flat with no header and have all the values for row 1, then row 2, etc
% the images and filters therefore need to be transposed for viewing in Matlab
% FilterConv3 replaces conv2() with fftolamopt2() from http://www.mathworks.com/matlabcentral/fileexchange/12171-ffw-fastest-filtering-in-the-west/all_files
% FilterConv4 replaces FilterConv3 by adding a single output file with the extension .filt for all the filtered images.

function [] =  FilterConv4(image_name, base_dir, Triples)
 
% close all hidden;
% format compact;

FileTypes = 2; % 1 for separate files for each of the filtered outputs (used for VisNetMat and VisNet C)
               % 2 for a single .filt file for MODE 3 of VisNetM3 and for VisNetMat1C.   3 for both file types

Display = 0; % set to 1 to display the images involved in the filtering; 0 for no display
% if >0, be prepared to touch any key to move on from the pause
% set to 2 to display also the neg and pos filtered images
% set to 3 to display asymmerically filtered images only valid with Triples = 1
%Triples = 0; % set to 1 to set the filtered images to have symmetric, asymmetric-low, asymmetric_high values in each row
% if Triples = 0, only the symmetric values will be output
%base_dir = '/home/erolls/VisNetUtilities/gaborL2/' % used for development
%base_dir = '/home/erolls/imagedata/imagesL3/'
%base_dir = '/home/twebb/src/imagedata/'
image_dir = strcat(base_dir, image_name, '.filtered/')
%filenameBase = '/home/erolls/imagedata/gaborL3/gabor'; % this is for the gabor filters produced by gaborfast512.c, e.g. gabor
%filenameBase = 'C:\VisNet\VisNetUtilities\gaborL3test\gabor2'; % this is for the gabor filters produced by gaborfast512.c, e.g. gabor
filenameBase = 'gabor'; % if the Gabor filters are in the current directory

mkdir(image_dir)

fig = 10;

Fsize = 512; % the padded size of the image, and the size of the filtered images
Imsize = Fsize/2; % the size of the input image
recon = zeros(Fsize, Fsize); % we will add all the symm filtered images here to see how well it reconstructs the original image, as a data check
recona = zeros(Fsize, Fsize); % we will add all the asymm filtered images here to see how well it reconstructs the original image, as a data check

% open and read in the image to be filtered
f = fopen(strcat(base_dir, image_name), 'r'); % a
tmp = fread(f, 'uint8=>uint8');
fclose(f);
% save the image to the .filtered directory with the added extension .original256
f_orig = fopen(strcat(image_dir, image_name,'.original256'), 'w'); 
tmp8 = uint8(tmp);
fwrite(f_orig, tmp8);
fclose(f_orig);

if FileTypes == 2 || FileTypes == 3
    % open the single file in the .filtered directory to contain all the filtered data with extension .filt
    f_filt = fopen(strcat(image_dir, image_name,'.filt'), 'w'); 
end

Image = reshape(single(tmp), [Imsize Imsize]); % for 512x512=262144 bytes
%Image = reshape(tmp, [256 256]); % for 256x256=65536 bytes
ImageMean = mean2(Image);
ImageStd = std2(Image);
fprintf('Image Mean=%3.1f  Std=%3.1f\n', ImageMean, ImageStd);
if (Display >= 1)
    fig = fig + 1;
    figure(fig);
    imagesc(Image'); title('Original Image');
    colormap(gray);
    pause;
% test only: add a black row below the 127 mean level
%    for (col=10: 246)
%        Image(220,col) = 0;
%        Image(221,col) = 0;
%    end
%    imagesc(Image); title('Original Image with 0 row');
%    colormap(gray);
%    pause;
end

% Make a mask to smooth the edges of the inserted image
SmoothMask = zeros(Imsize, Imsize);
SW = 10;
hpi = 2.0 * atan(1.0);
for row = 1:Imsize
    for col = 1:Imsize
        SmoothMask(row, col) = 1;
        x = 0;
        y = 0;
        if (row < SW) 
            x = SW - row;
        end
        if (row >= Imsize - SW)
            x = row + SW - Imsize + 1;
        end
        if (col < SW)
            y = SW - col;
        end
        if (col >= Imsize - SW)
            y = col + SW - Imsize + 1;
        end
        if ((x > 0) || (y > 0))
            r = sqrt(x * x + y * y);
            if (r > SW)
               r = SW;
            end
            SmoothMask(row, col) = 0.5 + 0.5 * sin(hpi * (SW - 2 * r) / SW);
        end
    end
end

% Now smooth the edges of the image to 0 so that there are no sharp edges in the paddedImage
paddingColor = 127; % Default is 127, but for ALOI set to 0
Image2 = Image - paddingColor;
Image2 = Image2 .* SmoothMask;
Image2 = Image2 + paddingColor;
if (Display >= 1)
    fig = fig + 1;
    figure(fig);
    imagesc(SmoothMask'); title('Smooth Mask');
    colormap(gray);
    fig = fig + 1;
    figure(fig);
    imagesc(Image2'); title('Smoothed Image');
    colormap(gray);
    MeanOfImage = mean2(Image2)
    pause;
end



% Create a new image of twice the dimensions
% and copy the original image into the center.
% Make padded image, set to provided background color
paddedImage = paddingColor * ones(Fsize,Fsize);
top = Fsize/2 - Imsize/2;
left = Fsize/2 - Imsize/2;
paddedImage(left + (1:Imsize), top + (1:Imsize)) = Image2;

if (Display >= 1)
    fig = fig + 1;
    figure(fig);
    imagesc(paddedImage'); title('Padded Image');
    meanVal = mean2(paddedImage)
    stdVal = std2(paddedImage)
    maxVal = max(max(paddedImage))
    minVal = min(min(paddedImage))
    colormap(gray);
    pause;
end


% produce the 128x128 padded image used by VisNet for its image display. This assumes Image is 256x256
%sub_sampled_image = subsampling(paddedImage, 4);
sub_sampled_image = imresize(paddedImage, 0.25, 'bicubic');

if (Display >= 1)
    fig = fig + 1;
    figure(fig);
    imagesc(sub_sampled_image'); title('Subsampled Padded Image 128x128');
    colormap(gray);
    pause;
end
f_sub = fopen(strcat(image_dir, image_name,'.image'), 'w');
sub_sampled_image = uint8(reshape(sub_sampled_image,[1 128 * 128]));
fwrite(f_sub, sub_sampled_image);
fclose(f_sub);

paddedImage = paddedImage - mean2(paddedImage);
if (Display >= 1)
    fig = fig + 1;
    figure(fig);
    imagesc(paddedImage'); title('Padded Image mean 0');
    meanVal = mean2(paddedImage)
    stdVal = std2(paddedImage)
    maxVal = max(max(paddedImage))
    minVal = min(min(paddedImage))
    colormap(gray);
    pause;
end

% For fftw library
load fftexecutiontimes
%t = cputime;
tic

%filenameBase = strcat(base_dir, 'gabor/gabor'); % this is for the gabor filters produced by gaborfast512.c, e.g. gabor

% List of frequencies and angles to iterate over
freqList = [1 2 4 8]; % from 28 Jan 2013
% freqList = [1 2 4 8 16 32 64 128];
%freqList = [128];
%freqList = [1 2 4 8]; % used for 32x32 VisNet
symAngList = [ 0 45 90 135]; % the symmetrical filters
asymAngList = [ 0 45 90 135 180 225 270 315]; % the asymmetrical filters

% load each Gabor filter and convolve each with the image
index = 1;
if (Display >= 1)
    fig = fig + 1;
    figure(fig);
end
for freq = freqList
    % Optimized parameters for 2D convolution
    %opt = detbestlength2(FFTrv,FFTiv,IFFTiv,size(paddedImage),size(sym_filterarray),isreal(paddedImage),isreal(sym_filterarray));

    for ang = symAngList
        %freq
        %ang
        % now open the next symmetric filter and do the convolution
        strcat(filenameBase, '.', int2str(freq), '.', int2str(ang),'.dat')
        fid = fopen( strcat(filenameBase, '.', int2str(freq), '.', int2str(ang),'.dat') );
        sym_filterarray = fread(fid, [8*freq*8*freq], 'single'); 
        fclose(fid);
        if (Display >= 1)
            Filter_sym_meanVal = mean(sym_filterarray)
            sym_stdVal = std(sym_filterarray)
            sym_maxVal = max(sym_filterarray)
            sym_minVal = min(sym_filterarray)
            % size(sym_filterarray)
        end
        sym_filterarray = reshape(sym_filterarray, freq*8, freq*8);
        %size(sym_filterarray)
        %sym_tmp = conv2(paddedImage, sym_filterarray, 'same');
        %sym_tmp = convolve2(paddedImage, sym_filterarray, 'same', 0.5);
        opt = detbestlength2(FFTrv,FFTiv,IFFTiv,size(paddedImage),size(sym_filterarray),isreal(paddedImage),isreal(sym_filterarray));
        sym_tmp = fftolamopt2(paddedImage, sym_filterarray, opt, 'same');
        sym_tmp = sym_tmp / freq; % scale filtered output by frequency to keep it in range -255 to 255
        recon = recon + sym_tmp;

        if (Triples == 1) % output symmetric, asymmetric_low, asymmetric_high triples in each row

            % now do the same for the two corresponding asymmetric filters
            fid = fopen( strcat(filenameBase, 'a.', int2str(freq), '.', int2str(ang),'.dat') );
            asym_low_filterarray = fread(fid, [8*freq*8*freq], 'single'); 
            fclose(fid);
            if (Display >= 1)
                Filter_asym_low_meanVal = mean(asym_low_filterarray)
                asym_low_stdVal = std(asym_low_filterarray)
                asym_low_maxVal = max(asym_low_filterarray)
                asym_low_minVal = min(asym_low_filterarray)
                % size(asym_low_filterarray)
            end
            asym_low_filterarray = reshape(asym_low_filterarray, freq*8, freq*8); % 30 Jun 2021
            % size(asym_low_filterarray)
            %asym_low_tmp = conv2(paddedImage, asym_low_filterarray, 'same');
            %asym_low_tmp = convolve2(paddedImage, asym_low_filterarray, 'same', 0.5);
            asym_low_tmp = fftolamopt2(paddedImage, asym_low_filterarray, opt, 'same');
            asym_low_tmp = asym_low_tmp / freq;
            recona = recona + asym_low_tmp;% comment out to exclude half the asymm filters from the reconstruction        
            
            fid = fopen( strcat(filenameBase, 'a.', int2str(freq), '.', int2str(ang+180),'.dat') );
            asym_high_filterarray = fread(fid, [8*freq*8*freq], 'single'); 
            fclose(fid);
            if (Display >= 1)
                Filter_asym_high_meanVal = mean(asym_high_filterarray)
                asym_high_stdVal = std(asym_high_filterarray)
                asym_high_maxVal = max(asym_high_filterarray)
                asym_high_minVal = min(asym_high_filterarray)
                % size(asym_high_filterarray)
            end
            asym_high_filterarray = reshape(asym_high_filterarray, freq*8, freq*8); % 30 June 2021
            % size(asym_high_filterarray)
            %asym_high_tmp = conv2(paddedImage, asym_high_filterarray, 'same');
            %asym_high_tmp = convolve2(paddedImage, asym_high_filterarray, 'same', 0.5);
            asym_high_tmp = fftolamopt2(paddedImage, asym_high_filterarray, opt, 'same');
            asym_high_tmp = asym_high_tmp / freq;
            recona = recona + asym_high_tmp; % comment out to exclude half the asymm filters from the reconstruction        
        
        end
        % Now create separate arrays for the positive, and the abs of the negative values
        % for a grayscale image 0,255, the filtered output will now be -127,127
        % set up the separate arrays for p and n, and clip at 255
        % some scaling will be necessary, especially for the low frequencies
        sym_tmp_pos = sym_tmp;
        sym_tmp_pos(sym_tmp < 0) = 0;
        sym_tmp_pos(sym_tmp > 255) = 255;

        sym_tmp_neg = sym_tmp;
        sym_tmp_neg(sym_tmp > 0) = 0;
        sym_tmp_neg = abs(sym_tmp_neg);
        sym_tmp_neg(sym_tmp < -255) = 255;
        
        if (Triples == 1) % output symmetric, asymmetric_low, asymmetric_high triples in each row
            % Asymmetric Low
            asym_low_tmp_pos = asym_low_tmp;
            asym_low_tmp_pos(asym_low_tmp < 0) = 0;
            asym_low_tmp_pos(asym_low_tmp > 255) = 255;

            asym_low_tmp_neg = asym_low_tmp;
            asym_low_tmp_neg(asym_low_tmp > 0) = 0;
            asym_low_tmp_neg = abs(asym_low_tmp_neg);
            asym_low_tmp_neg(asym_low_tmp < -255) = 255;

            % Asymmetric high
            asym_high_tmp_pos = asym_high_tmp;
            asym_high_tmp_pos(asym_high_tmp < 0) = 0;
            asym_high_tmp_pos(asym_high_tmp > 255) = 255;

            asym_high_tmp_neg = asym_high_tmp;
            asym_high_tmp_neg(asym_high_tmp > 0) = 0;
            asym_high_tmp_neg = abs(asym_high_tmp_neg);
            asym_high_tmp_neg(asym_high_tmp < -255) = 255;
        end
        
        output_pos = zeros(Fsize,Fsize); % for the filtered output positive values
        output_neg = zeros(Fsize,Fsize); % for the filtered output abs of neg values
        if (Triples == 0) % output symmetric values only
            output_pos = sym_tmp_pos;
            output_neg = sym_tmp_neg;
        end
        if (Triples == 1) % output symmetric, asymmetric_low, asymmetric_high triples in each row
            for row = 1:Fsize
                for col = 1:Fsize
                    if (mod(col, 3) == 0)
                        output_pos(row, col) = sym_tmp_pos(row,col);
                        output_neg(row, col) = sym_tmp_neg(row,col);
                    end
                    if (mod(col, 3) == 1)
                        output_pos(row, col) = asym_low_tmp_pos(row,col);
                        output_neg(row, col) = asym_low_tmp_neg(row,col);
                    end
                    if (mod(col, 3) == 2)
                        output_pos(row, col) = asym_high_tmp_pos(row,col);
                        output_neg(row, col) = asym_high_tmp_neg(row,col);
                    end
                end
            end
        end
        
        output_pos = uint8(reshape(output_pos, 1, Fsize*Fsize));
        output_neg = uint8(reshape(output_neg, 1, Fsize*Fsize));
        
        if FileTypes == 1 || FileTypes == 3
            fid_pos = fopen (strcat(image_dir,image_name, '.', int2str(freq), '.', int2str(ang),'.p'),'w');
            fwrite(fid_pos, output_pos);
            fclose(fid_pos);
            fid_neg = fopen (strcat(image_dir,image_name, '.', int2str(freq), '.', int2str(ang),'.n'),'w');
            fwrite(fid_neg, output_neg);
            fclose(fid_neg);
        end
        if FileTypes == 2 || FileTypes == 3
            fwrite(f_filt, output_pos);
            fwrite(f_filt, output_neg);
        end
        


        
        if (Display == 1)
            freq;
            ang;
            % size(sym_tmp)
            subplot(length(freqList), length(symAngList), index);
            set(gca,'xtick',[],'ytick',[]);
            imagesc(sym_tmp');  
            meanVal = mean2(sym_tmp)
            stdVal = std2(sym_tmp)
            maxVal = floor(max(max(sym_tmp)))
            minVal = floor(min(min(sym_tmp)))
            title({['Ang' num2str(ang) ' Max ' num2str(maxVal) ' Min ' num2str(minVal)  ]});
        end

        if (Display == 3) % show the asymm filter activity
            freq;
            ang;
            subplot(length(freqList), length(symAngList)*2, index);
            set(gca,'xtick',[],'ytick',[]);
            imagesc(asym_low_tmp');  
            meanVal = mean2(asym_low_tmp)
            stdVal = std2(asym_low_tmp)
            maxVal = floor(max(max(asym_low_tmp)))
            minVal = floor(min(min(asym_low_tmp)))
            title({['Ang' num2str(ang) ' Max ' num2str(maxVal) ' Min ' num2str(minVal)  ]});
            
            subplot(length(freqList), length(symAngList)*2, index+1);
            set(gca,'xtick',[],'ytick',[]);
            imagesc(asym_high_tmp');  
            meanVal = mean2(asym_high_tmp)
            stdVal = std2(asym_high_tmp)
            maxVal = floor(max(max(asym_high_tmp)))
            minVal = floor(min(min(asym_high_tmp)))
            title({['Ang' num2str(ang+180) ' Max ' num2str(maxVal) ' Min ' num2str(minVal)  ]});
            index = index + 1;
        end

        
% display firing rates of the filtered outputs, which provide the inputs to VisNet
% Angles are arranged horizontally, frequencies vertically with high at the top
% For each angle, there is an image of the firing in array where the filter output was positive,
% and then immediately to the right for the same angle where the filter output was negative 
        if (Display == 2)
            subplot(length(freqList), length(symAngList)*2, index);
            set(gca,'xtick',[],'ytick',[]);
            imagesc(sym_tmp_pos');  
            colormap(gray);
            meanValpos = mean2(sym_tmp_pos)
            stdValpos = std2(sym_tmp_pos)
            maxValpos = floor(max(max(sym_tmp_pos)))
            minValpos = min(min(sym_tmp_pos))
            title({['Ang' num2str(ang) ' Pos Max ' num2str(maxValpos) ' Min ' num2str(minValpos)  ]});
            %pause

            subplot(length(freqList), length(symAngList)*2, index+1);
            set(gca,'xtick',[],'ytick',[]);
            imagesc(sym_tmp_neg'); 
            colormap(gray);
            meanValneg = mean2(sym_tmp_neg)
            stdValneg = std2(sym_tmp_neg)
            maxValneg = floor(max(max(sym_tmp_neg)))
            minValneg = min(min(sym_tmp_neg))
            title({['Ang' num2str(ang) ' Neg Max ' num2str(maxValneg) ' Min ' num2str(minValneg)  ]});
            index = index + 1;
        end
        
    index = index + 1;
    end
end

if FileTypes == 2 || FileTypes == 3        
    fclose(f_filt); % the single output file
end
% colormap(gray);

%t = cputime - t;
%disp(['Time for convolutions: ' num2str(t)]);
toc

if (Display >= 1)
    pause;
    % display the gabor filters
    symFilters = cell(1, length(freqList)*length(symAngList)); % used just for the display
    %asymFilters = cell(1, length(freq)*length(asymAngList));
    asymFilters = cell(1, length(freqList)*length(asymAngList));
    fig = fig + 1;
    figure(fig); % display the symmetrical filters
    %title('Symmetric filters'); hold on;
    colormap default;
    index = 1;
    for freq = freqList
        for ang = symAngList
            fid = fopen( strcat(filenameBase, '.', int2str(freq), '.', int2str(ang),'.dat') );
            symFilters{index} = fread(fid, [8*freq 8*freq], 'single'); 
            fclose(fid);
            
            subplot(length(freqList), length(symAngList), index);
            imagesc(symFilters{index}');
            %image(symFilters{index});
            set(gca,'xtick',[],'ytick',[]);
            index = index + 1;
            sgtitle('Symmetric filters')
        end
    end
    pause

      
    if (Triples == 1)
        fig = fig + 1; % display the asymmetrical filters
        figure(fig);
        colormap default;
        % title('Asymmetric filters'); hold on;
        index = 1;
        for freq = freqList
            for ang = asymAngList
                fid = fopen( strcat(filenameBase, 'a.', int2str(freq), '.', int2str(ang),'.dat') );
                asymFilters{index} = fread(fid, [8*freq 8*freq], 'single'); 
                fclose(fid);
                subplot(length(freqList), length(asymAngList), index);
                imagesc(asymFilters{index}');
                set(gca,'xtick',[],'ytick',[]);
                index = index + 1;
                sgtitle('Asymmetric filters')
            end
        end
        %title('Asymmetric filters'); hold on;
    end
end

if (Display >= 1)
    fig = fig + 1;
    figure(fig);
    imagesc(recon'); title('Reconstructed Image from symmetrical filtered images');
    colormap(gray);
        
    if (Triples == 1)
        fig = fig + 1;
        figure(fig);
        imagesc(recona'); title('Reconstructed Image from asymmetrical filtered images');
        colormap(gray);

        fig = fig + 1;
        figure(fig);
        reconboth = recon + recona;
        imagesc(reconboth'); title('Reconstructed Image from both symm and asymm filtered images');
        colormap(gray);
    end
end

    
    