%% ========================================================================
% Demo codes for Testing for Super Resolution
% =========================================================================
% Code Developed by Pulak Purkait 
% ------------------------------
% ECSU, Indian Statistical Institute Kolkata
% -------------------------------------------
% For any questions/bug related issues, send email to: pulak_r@isical.ac.in
% Referece 
%      P. Purkait and B chanda, Fuzzy-Rule Based Approach for Single Frame 
%                                           Super Resolution, FuzzIEEE 2013
% =========================================================================

% Lim = imread('001_input.png');  % Input image 
clear;
Lim = imread('input_003.png');  % Input image 

rs_fact = 4;                    % resolution Factor

chnl = size(Lim, 3);
if chnl == 3,
    im_ycbcr = rgb2ycbcr(Lim);
    lImL = double(im_ycbcr(:, :, 1));
else
    lImL = double(Lim);
end
lImO = lImL;

%% Load Dictionaries and the parametres trained for Super Resolution  
dict_path = 'Dictionary/D_75_7_1.6.mat';
load(dict_path);

D(isnan(D)) = 0;

%% Parameters for EM algorithm

% par.maxit  = 25; % maximum number of iterations of EM
% par.Leps   = 0.1; % EM stopping condition: minimum change in log-likelihood
par.doplot = 0;   % Ploting verbosity

par.ovrlap = 2;     % Document overlap size 
par.proj_mxItr = 10;    % Number of iteration for Back-Projection 

for ind = 1:round(log(rs_fact)/log(par.upscale))

    % Upscale the image 
    
    lIm = imresize(lImO, par.upscale, 'bicubic');   % Image initialize by Bicubic interpolation 
    [nrow, ncol] = size(lIm);

    nrow2=ceil((nrow - par.ovrlap)/(par.patch_size-par.ovrlap))*(par.patch_size-par.ovrlap)+par.ovrlap;
    ncol2=ceil((ncol - par.ovrlap)/(par.patch_size-par.ovrlap))*(par.patch_size-par.ovrlap)+par.ovrlap;

    lIm2 = lIm([1:nrow, nrow*ones(1, nrow2-nrow)], [1:ncol, ncol*ones(1, ncol2-ncol)]);

    [WrdsA, Xvar, Mns, Xorg] = test_sample_patches(lIm2, par.patch_size, par.ovrlap);

    % Selective Patch Processing 
    Wrds = WrdsA(1:end-par.patch_size^2, :); 
    Hptchs = WrdsA(end-par.patch_size^2+1:end, :);
    idx = sum(Wrds.^2) > 20000; 
    Wrds = Wrds(:, idx); 
    Wrds = Wrds./repmat(sqrt(sum(Wrds.^2))+eps, size(Wrds, 1), 1); 

    % % Calculate topic for each visual words
    tic
    A = computeA(Wrds', NewC, NewSig);
    Hptchs(:, idx) = (A*D)'; 
    toc
    Hptchs = Hptchs./repmat(sqrt(sum(Hptchs.^2))+eps, size(Hptchs, 1), 1); 
    
    OutVar = 1.2*Xvar; % Output patch varience.
    Hptchs = Hptchs.*repmat(OutVar,[size(Hptchs,1) 1]) + repmat(Mns, [size(Hptchs,1) 1]);

    % % Generating HR Image

    hIm = zeros([nrow2, ncol2]);
    hIm2 = hIm; 

    count = zeros([nrow2, ncol2]);
    x = 1:par.patch_size-par.ovrlap:nrow2-par.patch_size+1;
    y = 1:par.patch_size-par.ovrlap:ncol2-par.patch_size+1;

    [XX,XY] = meshgrid(x,y);

    xrow = XX(:);
    ycol = XY(:);

    % compute the The high resolution patches using topic model

    cntn = 0;

    for ii = 1:length(xrow)

        row = xrow(ii)-1;
        col = ycol(ii)-1;

        hIm( row+1: (row+par.patch_size), col+1: (col+par.patch_size)) = ...
            hIm( row+1: (row+par.patch_size), col+1: (col+par.patch_size)) ...
            + reshape(Hptchs(:, ii), [par.patch_size, par.patch_size]);

        count( row+1: (row+par.patch_size), col+1: (col+par.patch_size)) 	= ...
            count( row+1: (row+par.patch_size), col+1: (col+par.patch_size)) + ones(par.patch_size);

    end

    Id = count>0;
    hIm(Id) = hIm(Id)./count(Id);
    hImm = hIm(1:nrow, 1:ncol);

    if ind ~= 3
        hImm = backprojection(hImm, lImL, par.proj_mxItr);
    end

    lImO = hImm; % Initialize LR image 
end

if chnl == 3
    Lim = imresize(Lim, size(hImm), 'bicubic');
    im_ycbcr = rgb2ycbcr(Lim);
    im_ycbcr(:, :, 1) = hImm;
    hImm = ycbcr2rgb(im_ycbcr);

else
    Lim1 = imresize(Lim, size(hImm), 'bicubic');
    imshow(uint8(hImm));
    figure, imshow(uint8(Lim1));
%     figure, imshow(uint8(255*clrplot));
end

imwrite(uint8(hImm), 'FuzzyRulebase.png');
imwrite(uint8(Lim), 'Bicubic.png');


