
function fix_images()

close all
clear all

% nekde taky udelat resuduum povudni minus opravený !!!

global W H Gap path

path = 'photo_side/'  % load/save path 

% images parameters (constants)
W = 96 ;
H = 336;
Gap = 16;



clear CorrectedRGB 
for i = 1:3
    img_all =load_img();
    data =  load_data(img_all, i, true);
%     keyboard
    range = 1:size(data,2); % process only some images 
    data = data(:,range);
    OrigRGB(:,:,i) =  double(reshape(data, H, []));
    CorrectedRGB(:,:,i) =  main(img_all, i, range);
%     imagesc(CorrectedRGB(:,:,i))
    
end

keyboard



% parameters used becuase of the nonlinearity the photos
a = 43000;%+-600;
b = 140;%+-3;
c = 51/1460;

OrigRGB         = OrigRGB.*a./(1/c+b*OrigRGB);
CorrectedRGB    = CorrectedRGB.*a./(1/c+b*CorrectedRGB);
OrigRGB         = uint8(OrigRGB);
CorrectedRGB    = uint8(CorrectedRGB);

% figure
imshow(CorrectedRGB)
 
N = size(CorrectedRGB,2);
imwrite(CorrectedRGB,[path 'CorrectedRGB.png']);
imwrite(OrigRGB,[path, 'OrigRGB.png']);

try
    diff = double(CorrectedRGB) - double(OrigRGB(:,1:N,:));
    diff = diff(:,:,1);
    diff(:,sum(sum(OrigRGB(:,1:N,:),3)) == 0) = nan;
    imagesc(diff(:,:))
    figure
    imagesc(diff(:,7940:9070)); colorbar;   
    axis image
    axis off
    set(gcf, 'PaperType', 'A1')

%     imwrite(diff,[path, 'diff.png']);
end

keyboard

  
end


function data_new = main(img_all, color, range)

    % load data for one color and call all filters and gap filling
    

    data =  load_data(img_all,  color, false);  % barva !!
    data =  cutoff_data(data);

    N_img =  size(data,2);

%     keyboard
    
%%%%%%%%%  REMOVE REFLECTIONS %%%%%%%%%%%%%%%%%%%%%%
    % use the strongest filter only for the 2 first vectors
    % substract them from the data 
    for i = 1:2; [~, vec{i}] = total_filter(data, i ); end
    for i = 1:2; s{i} = ((data'*vec{i}')*vec{i})'; end
    for i = 1:2; data = data - s{i}; end

    
    concat_images(data, range, true);  % show results
    data =  cutoff_data(data);  % remove simple reflections
    [data, ~] = total_filter(data, 1 );  % strong smoothing
    N_filt = min(10, ceil(N_img/10));
    data = poly_filter(data, 1:N_filt, 3);  % minor smoothing
    concat_images(data, range, true);

    % add the data back
    for i = 1:2; data = data + s{i}; end

    % try to remove it once more
    data =  cutoff_data(data); % remove simple reflections
    data =  total_filter(data,1); % minor smoothing ( once more)
   
%%%%%%%%%%%%%%  FILL GAPS %%%%%%%%%%%%%%%

%     data_new = concat_images(data, range);  % fill gaps with nans
    data_new = gap_filling(data, range); % fill shots with "data"

end



function img_all =load_img()
    % load data from harddisk

    global  path
    name = [path, 'cache.mat'];
    try
        i = load(name);
        img_all = i.img_all;
    catch
        disp('loading ...')
        k = 1;
        for j = 7000:9000  % load all images from the folders ....
            for i =   [2000:3000, 8000:9000]
                p = [path, sprintf('%i',j), '-0000',  sprintf('%04.0f',i ) ,'.png'];
                if exist(p, 'file')
                    img = imread(p);
                    if ndims(img) == 3 
                        img_all{k} = img;
                        k=k+1
                    end
                end
            end
        end
        disp('cache saved  ..')
        save(name, 'img_all')
    end
end

function  data =  load_data(img_all, color, add_gaps)
    % process the loaded images

    global W H N_img Gap 

    N_img = length(img_all);

    for i = 1:N_img   
        img = double(img_all{i});
        %  convert the intensity to the linear scale
        a = 43000;%+-600;
        b = 140;%+-3;
        c = 51/1460;
        img = img./(a-b*img)/c; %normované od 0 do 1

        if ~add_gaps  % add black gaps to the images
            image(:,:,i) = sum(img(:,:,color),3);
        else
            image(:,:,i) = [sum(img(:,:,color),3), zeros(H,Gap)];
        end
    end

    data =  reshape(image, (W+add_gaps*Gap)*H,[])+max(image(:)*1e-6);  % make all pixels nonzero

end


function img_concat = concat_images(data, range, plot)
    % connect images without gap filling 

    global W H Gap
    
   N_img = length(range);
   for i = 1:N_img
        img_concat(:,(i-1)*(W+Gap) + (1:W)) = reshape(data(:,range(i)), H,W);
   end
   
   if nargin > 2
     figure;  imagesc(img_concat)
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function  data =  cutoff_data(data)
    global W H 

    % remove reflections from the top of the images
    data = reshape(data, H,W,[]);
    data(1:25,:) = 0;
    data = reshape(data, H*W,[]);
    data(data < 0) = 0;   
end

function [data, x] = total_filter(data, N_vec)

    % remove all time resolution from the N_vec-th vector

    global W H % decompose the images, each image is one vector
    [u,s,v] = svd(data, 'econ');
    base = reshape(u, H,W,[]);
    filtered = mean(base(:,30:70,N_vec),2)'; % use only centre of the image 
    filtered = smooth(filtered,20);  % smooth in the "y" direction
%     figure
%     imagesc(filtered)
%     keyboard
    x = repmat(filtered, W,1)';

    u(:,N_vec) = x(:);  % replace is the original vector
    data = u*s*v';   % composition of the original matrix

end


function data = poly_filter(data, ind_vec, poly)
    global W H 

    % strong smoothing in order to remove static artifacts
    
    [u,s,v] = svd(data, 'econ');  % decompose the images, each image is one vector

    base = reshape(u, H,W,[]);

    for i = ind_vec
        for j = 1:H
            % brutal smoothing of the vectors
            p = polyfit(1:W, base(j,:,i), poly);
            d  = polyval(p,1:W);
            base(j,:,i) = d(:);
        end
    end
    u = reshape(base, W*H, []);
    data = u*s*v';
    
end



function data_new = gap_filling(data_orig, range)
    % fill the gaps with a good estimation 

    global W H Gap

    N_pic = length(range);
        
    % reshape the images to the original size of the photos 
    for i = 1:N_pic;
        img_all_conc(:,(i-1)*W + (1:W)) = reshape(data_orig(:,range(i)), H,W);
    end
% keyboard
    % decompose the images (each image column is one vector)
     
    img_all_conc = img_all_conc + 0.01;  % important fix to avoid problems with SVD !!!!
    
    [u,s,v] = svd(img_all_conc, 'econ');
 
%      imagesc(u*s*v')
 
   
    time = v(1:N_pic*W,:);


    N=H*2+Gap;
    Nfft=2^nextpow2(N); % speed up the fft

    % add gaps 
    time_new = nan(Nfft, H);
    for i = 1:N_pic
        time_new((i-1)*(W+Gap)  + (1:W) , :) = time((i-1)*W  + (1:W) , :);
    end

    % keyboard
    time_tmp = time_new;

    N=W+Gap;
    Nfft=2^nextpow2(N);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  FILL GAPS %%%%%%%%%%%%%%%%%%%%%%%
    for i = 1:H
        i
        tic
        for j = 1:N_pic-1
            data = nan(Nfft, 1);
            data(1:N) = [time((j-1)*W  + (W/2:W) , i); nan(Gap,1) ; time(j*W  + (1:W/2-1) , i)];
%             data_0 = data;
            if i < 50  % only the most important vectors 

                ind = find(isnan(data));
                for k = [-1,1]
                    ind = unique([ind+k, ind]);
                end
                ind = ind(ind > 0 & ind < Nfft);
                data(ind) = nan;
                data_tmp = ifft(edft(data), Nfft);      % !!! the main woodoo !!! , quite slow ..
                if any(~isreal(data_tmp))
                    keyboard
                end
            else
%                 try
                data_tmp = interp1(find(~isnan(data)),data(~isnan(data)), 1:(W+Gap) );  % fast cheating for the less important vectors 
%                 catch
%                     keyboard
%                 end
            end
            time_new((j-1)*(W+Gap)  + W+ (1:Gap), i) = data_tmp(W/2+(1:Gap));
        end
        toc
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% keyboard


    data_new = u*s*time_new';

end
