% Mahmood Amintoosi, IUST 2009
% This is the main file for mapping and fusing
function [converged,convergedIter,rms_error,ssim_error,times] = run_homo_MAT_10_LM(n_iters,alg,snr,lrImName)

% NOTE: Since we used ~exist, if a variable not exist it will be set only
% in the first run.
if ~exist('n_iters')
    n_iters = 15;					% Number of gradient descent iterations
end
if ~exist('alg')
    alg.list = {'homo_fa_MAT' 'homo_fa_MAT_SSIM' 'homo_fa_lm_MAT_SSIM'};
    alg.name = {    'LK'       'LK-SSIM' 'LK-SSIM-LM'};
end
if ~exist('snr')
    snr = 20; %ICEE, snr=90
end
if ~exist('lrImName')
    lrImName = '0000'; %ICEE, snr=90
end

im1FilteringSW = 0; %0:No Filtering,  1: Low pass filtering of the first image
im2FilteringSW = 0; %0:No Filtering,  1: Low pass filtering of the first image
addNoiseToIm2SW = 1;
croppingIm1SW = 0;
freemanSW = 0;
showSW = 0;%3;
close all;
%% Reading image and template
% im1FileName='D:\PhDCourses\SuperResolution\Video_SR\HR_Images\Shanghai_Cropped7.jpg';
% im2FileName='D:\PhDCourses\SuperResolution\Video_SR\HR_Images\Shanghai_1.jpg';
% % im2FileName='Shanghai_qlty_70.jpg';
% desc = 'Shanghai';masir = desc; hrScale=.5; lrScale=2; im1FilteringSW = 0;

% im1FileName='E:\Data\Bistoon\ScannedPhotoes\Darius_scan0082.JPG';
% im2FileName='E:\Data\Bistoon\cam\Katibeh_Darius_half.JPG';%hrScale=2; lrScale=6;
% desc = 'Darius';masir = 'Darius_02'; hrScale=.5; lrScale=1.5;
% im1FilteringSW=1;croppingIm1SW = 0;

% im1FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\1b8\YaFatemeh_1b8_LR_4.tif';
% im2FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\YaFatemeh_LR_1.tif';
% desc = 'YaFatemeh';masir = desc; hrScale=2; lrScale=1;im1FilteringSW = 0;

% im1FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\Breakfast_LR_1_cropped.tif';
% % im2FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\Breakfast_snr10_LR_4.tif';
% im2FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\Breakfast_LR_4.tif';
% desc = 'Breakfast';masir = [desc '_1minus']; hrScale=1; lrScale=1;
% im1FilteringSW = 0; croppingIm1SW = 0;  im2FilteringSW=0;

% im1FileName='D:\PhDCourses\SuperResolution\SR_Init\Park\Park1_01.jpg';
% im2FileName='D:\PhDCourses\SuperResolution\SR_Init\Park\Park1_02_LR.jpg';
% desc='Sorsoreh'; masir = 'Sorsoreh'; hrScale=.5; lrScale=2; im1FilteringSW=0;

% im1FileName='D:\PhDCourses\Image\Peter_Kovesi_MatlabFns\Robust\example\testImages\Balash1.JPG';
% im2FileName='D:\PhDCourses\Image\Peter_Kovesi_MatlabFns\Robust\example\testImages\Hadis_Balash332.jpg';
% desc = 'Balash'; masir = 'Balash'; hrScale=2; lrScale=2;

% % im1FileName='D:\PhDCourses\Image\Peter_Kovesi_MatlabFns\Robust\example\testImages\Closing_HR.jpg';
% im1FileName='D:\PhDCourses\SuperResolution\Video_SR\HR_Images\Closing_HR_Cropped2.jpg';
% im2FileName='D:\PhDCourses\Image\Peter_Kovesi_MatlabFns\Robust\example\testImages\Closing_LR_4.tif';
% desc='Closing_HR4x_LR4';masir = 'Closing';hrScale=1; lrScale=2; im1FilteringSW=0; im2FilteringSW=0;

% im1FileName='D:\PhDCourses\Image\Peter_Kovesi_MatlabFns\Robust\example\testImages\DSC02607.JPG';
% im2FileName='D:\PhDCourses\Image\Peter_Kovesi_MatlabFns\Robust\example\testImages\DSC02609.JPG';
% desc = 'Longhua'; masir = desc; hrScale=.2; lrScale=.2;

% im1FileName='C:\Documents and Settings\Mahmood\Desktop\China\DSC02641.JPG';
% im2FileName='C:\Documents and Settings\Mahmood\Desktop\China\DSC02640.JPG';
% desc = 'Yanan'; masir = desc; hrScale=.25; lrScale=.35;

% im1FileName='C:\Documents and Settings\Mahmood\Desktop\China\DSC02877.JPG';
% im2FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\InBus_LR.bmp';
% desc = 'InBus';masir = desc; hrScale=.25; lrScale=2;im1FilteringSW = 0;

% % im1FileName='C:\Documents and Settings\Mahmood\Desktop\China\DSC02645.JPG';
% im1FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\Bayram_HR_Cropped.jpg';
% im2FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\Bayram_LR_Cropped2.bmp';
% desc = 'Bayram';masir = desc; hrScale=.25; lrScale=2;im1FilteringSW = 0;

% % For IUST ICEE2009
% im1FileName='E:\Data\Bistoon\CD\Neyzeh_dar.JPG';%Daste_Kamandar.jpg';%Sare_Kamandar.jpg';%
% im2FileName='E:\Data\Bistoon\CD\Afsaran.jpg';
% % For IUST ICEE2009
% % desc = 'Kamandar';masir = desc; hrScale=1; lrScale=2;im1FilteringSW = 0;
% % For Lahijan
% desc = 'Kamandar';masir = 'KamandarLM'; hrScale=1; lrScale=2;im1FilteringSW = 0;
% 
% im1FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\YuyuanGarden_HR_Cropped.jpg';
% im2FileName='D:\PhDCourses\SuperResolution\Video_SR\testImages\YuyuanGarden_LR.bmp';
% desc = 'YuyuanGarden';masir = desc; hrScale=.25; lrScale=2;im1FilteringSW = 0;

% %For Video LM
im1FileName='E:\Data\My_SR_Data\Tokyo\THR_Images\0218.bmp';%
im2FileName=['E:\Data\My_SR_Data\Tokyo\LR_Images_half\' lrImName '.bmp'];%
% im2FileName='E:\Data\My_SR_Data\Tokyo\LR_Images_half\0010.bmp';
desc = 'Tokyo';masir = desc; hrScale=1; lrScale=2;im1FilteringSW = 0;

if exist(masir,'dir')~=7
    mkdir(masir);
end
desc = [masir '\' desc];

filtWinSize = 5;

im1 = imread(im1FileName);
im1 = imresize(im1,hrScale); 
if im1FilteringSW==1
    im1 = imfilter(im1,fspecial('average',filtWinSize),'replicate'); %Low pass filter
end
im1Color = im1;
im1 = im2double(im1);
im1 = rgb2gray(im1);
% im1(1:floor(58*size(im1,1)/100),:) = [];
% SNR=20;%:10:20;
% for snrNo=1:length(SNR)
%     snr = SNR(snrNo);
im2 = imread(im2FileName);
if addNoiseToIm2SW==1
    im2 = add_noise(im2,snr);
end
if im2FilteringSW==1
%     im2 = medfilt2(im2,[5 5]);
%     im2 = wiener2(im2,[5 5]);
    im2 = imfilter(im2,fspecial('average',filtWinSize),'replicate'); %Low pass filter
end
im2ColorNearest = imresize(im2,lrScale,'nearest');
im2Color = imresize(im2,lrScale,'bicubic');
im2 = im2Color;
im2 = im2double(im2);
im2 = rgb2gray(im2);


if croppingIm1SW==1 && snrNo==1
    figure(1); clf;
    imshow(im1);
    rect = floor(getrect);
    im1 = im1(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3));
end

options.desc = desc;
options.masir = masir;
options.thresh = .20;
options.distRatio = .8;
options.fileName1 = sprintf('%s_%dx%d',getFileName(im1FileName),size(im1,1),size(im1,2));
options.fileName2 = sprintf('%s_%dx%d_snr%d',getFileName(im2FileName),size(im2,1),size(im2,2),snr);
options.showSW = showSW;
% methodName = 'Corners';
methodName = 'SIFT';
options.mask = true(size(im1,1),size(im1,2));
[im3, mappingParam ] = MAT_homog(im1,im2,methodName,options);
H = mappingParam.H;

if showSW ~= 0
    figure(11); %subplot(121); imshow(options.mask); subplot(122); imshow(mapOutOptions.mask)
    imshow(im3);%imshow(im1Transformed);
%     figure(3); figure(4);
    
    figure(9);
    imshow(im2);
end
options.showSW = showSW;
[im1Transformed,im1OnIm2] = Im1Transform(im1Color,im2Color,H,options);
if showSW ~= 0
    figure(11); %subplot(121); imshow(options.mask); subplot(122); imshow(mapOutOptions.mask)
    imshow(im1OnIm2);%imshow(im1Transformed);
end
tmplt = im1;
img = im2; imgColor = im2Color;
H = H./H(9);

if showSW ~= 0
    fileName = sprintf('%s_img.jpg',desc);
    imwrite(imgColor,fileName);
    fileName = sprintf('%s_tmplt.jpg',desc);
    imwrite(im1Color,fileName);
end

im1OnIm2_Fustruction = im1OnIm2;
% imwrite(im1OnIm2,'img1OnIm2BeforePA.jpg');
% save testClosing.mat img tmplt H
% Test parameters
verbose = 0;					% Show fitting?
clear fit times;

for algNo=1:length(alg.list)
		string = ['tic; fit{algNo} = ', alg.list{algNo}, '(img,tmplt,H,n_iters,verbose); t = toc;'];
		eval(string);
        times(algNo)=t;
end
%% Convergence Speedup
min_rms = min([fit{1}.rms_error]);
first_rms = fit{1}(1).rms_error;

% Here instead of considering the min_rms over all methods, we consider the
% FA method best rms
% for algNo = 2:length(alg.list)
%     min_rms = min([min_rms, [fit{algNo}.rms_error]]);
% end
min_rms = min([fit{1}.rms_error]);

nintyPercent_rms = 85*(first_rms-min_rms)/100;
convergingIter= ones(1,length(alg.list))*(n_iters+1);
converged = zeros(1,length(alg.list));

% The corresponding element with the first iter number, which algorithm
% converged, is saved in convergedIter
convergedIter = false(length(alg.list),n_iters);

for algNo = 1:length(alg.list)
    %if last rms is greater than first, the algorithm diverged
    if fit{algNo}(end).rms_error>first_rms, continue, end
    tmp = [fit{algNo}.rms_error]<(first_rms-nintyPercent_rms);
    idx = find(tmp>0);
    if ~isempty(idx),    convergingIter(algNo) = idx(1); end
end

% n_iters is added for the case which all algorithms diverged
% converged indicates which algorithm converged faster.
minIter = min([n_iters,convergingIter]);
for algNo = 1:length(alg.list)
    converged(algNo) = convergingIter(algNo)<=minIter;
    if convergingIter(algNo)<=n_iters, convergedIter(algNo,convergingIter(algNo))=true; end
    rms_error{algNo} = [fit{algNo}.rms_error];
    ssim_error{algNo} = [fit{algNo}.ssim_error];
end
% converged

bar(times);
title('Run Time'); 
set(gca,'XTickLabel', alg.name);% ,'fontsize',12);
% Plotting errors
markers = {'rd--','bp-','ks--','g.-','+:','s:','d:'};
% RMS Error
h=figure(21); clf;
set(h,'DefaultLineLineWidth',1);
set(gcf,'name',sprintf('%s, snr=%d',desc,snr));
% subplot(211);
hold on;
for algNo=1:length(alg.list)
    plot([fit{algNo}.rms_error],markers{algNo});
end
xlabel('Iteration No.');
ylabel('RMS Error');
% title('Example RMS-Error'); 
legend(alg.name);
set(gca,'XTick', 0:n_iters);

% xtmp = 0:n_iters;
% ytmp = ones(1,length(xtmp)).*first_rms-nintyPercent_rms;
% h2=plot(xtmp,ytmp,'m:');
% set(h2,'LineWidth',2);
% hold off;
% print('-depsc','example_rms_error.eps');
print_EPS(h,sprintf('%s_rms_error_%2d.eps',desc,snr),[4.3 3]);

% SSIM Error
h=figure(22);clf;
set(h,'DefaultLineLineWidth',1);
hold on;
for algNo=1:length(alg.list)
    if ~isfield(fit{algNo},'ssim_error'), continue, end
    plot([fit{algNo}.ssim_error],markers{algNo});
end
xlabel('Iteration No.');
ylabel('SDIS Error = -MSSIM');
%title('Example SSIM-Error'); 
legend(alg.name);
hold off;
% print('-depsc','example_ssim_error.eps');
print_EPS(h,sprintf('%s_ssim_error_%2d.eps',desc,snr),[4.3 3]);
% %% Cropping Region 
% nBand = length(size(im2Color));
% figure(1);
% imshow(imgColor); title('The original LR image');
% imwrite(imgColor,sprintf('%s_LR.jpg',desc));
% rect = floor(getrect);
% temp = imgColor(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3),:);
% imwrite(temp,sprintf('%s_LR_cropped.jpg',desc));
% 
% imwrite(im2ColorNearest,sprintf('%s_LR_Nearest.jpg',desc));
% temp = im2ColorNearest(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3),:);
% imwrite(temp,sprintf('%s_LR_Nearest_cropped.jpg',desc));
% 
% if freemanSW==1
%     temp = imread(sprintf('%s_Freeman.jpg',desc));
%     freeman = imresize(temp,[size(imgColor,1),size(imgColor,2)],'bicubic');
%     temp = freeman(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3),:);
%     imwrite(temp,sprintf('%s_freeman_cropped.jpg',desc));
% end
% 
% im2ColorSyn = im2uint8(im1OnIm2_Fustruction);
% XFus = zeros(size(im2ColorSyn));
% for band=1:nBand
%     XFus(:,:,band)  = wfusimg(double(im2Color(:,:,band)),double(im2ColorSyn(:,:,band)),'db2',3,'img1','img2');
% end
% Fustruction = uint8(XFus);
% imwrite(Fustruction,sprintf('%s_Fustruction.jpg',desc));
% temp = Fustruction(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3),:);
% imwrite(temp,sprintf('%s_Fustruction_cropped.jpg',desc));
% %% 
% for algNo=1:length(alg.list)
% %     [val,idx]=min([fit{algNo}.rms_error]);
%     idx = round(times(1)/times(algNo)*n_iters); %For equal time
%     [im1Transformed,im1OnIm2] = im1Transform(im1Color,im2Color,fit{algNo}(idx).warp_p,options);
%     im2ColorSyn = im2uint8(im1OnIm2);
%     %Fusion
%     XFus = zeros(size(im2ColorSyn));
%     for band=1:nBand
%         XFus(:,:,band)  = wfusimg(double(im2Color(:,:,band)),double(im2ColorSyn(:,:,band)),'db2',3,'img1','img2');
%     end
%     XFus = uint8(XFus);
%     imwrite(XFus,sprintf('%s_XFus_%s.jpg',desc,alg.name{algNo}));
%     temp = XFus(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3),:);
%     imwrite(temp,sprintf('%s_XFus_%s_cropped.jpg',desc,alg.name{algNo}));
%     figure(algNo+1); % imshow(options.mask); subplot(122); imshow(mapOutOptions.mask)
% %     subplot(121);
% %     imshow(XFus);%im1OnIm2); 
%     title(sprintf('rms=%f , ssim=%f',fit{algNo}(idx).rms_error,fit{algNo}(idx).ssim_error));
%     
%     %Blending 
%     mask = rgb2gray(im1Transformed)~=0;
% %     %Filling the holes of mask
% %     mask = bwmorph(mask,'close',Inf); 
% %     %Reducing the size of mask for avoiding the borders of wfusimg
% %     se1 = strel('square',11); mask = imerode(mask,se1); 
% % 
% %     hrIm = im2double(rgb2gray(XFus));
% %     lrIm = im2double(im2);
% %     [hrIm,lrIm,mask,extended] = power_of_2_images(hrIm,lrIm,mask);
% %     se1 = strel('square',11);
% %     mask = imerode(mask,se1);
% %     [blendedIm,D] = BlendArbitrary(hrIm, lrIm, mask, 4);
% 
%     blendedIm = multiBandBlending(im1Transformed,im2Color,mask);
%     subplot(121);
%     imshow(blendedIm); title('Blending of transformed of HR');
%     imwrite(blendedIm,sprintf('%s_blend_%s.jpg',desc,alg.name{algNo}));
%     temp = blendedIm(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3),:);
%     imwrite(temp,sprintf('%s_blend_%s_cropped.jpg',desc,alg.name{algNo}));
%     
%     blendedIm = multiBandBlending(XFus,im2Color,mask);
%     subplot(122);
%     imshow(blendedIm); title('Blending of Fusion');
%     imwrite(blendedIm,sprintf('%s_blendXFus_%s.jpg',desc,alg.name{algNo}));
%     temp = blendedIm(rect(2):rect(2)+rect(4),rect(1):rect(1)+rect(3),:);
%     imwrite(temp,sprintf('%s_blendXFus_%s_cropped.jpg',desc,alg.name{algNo}));
% %     pause
% end
% 
% print('-djpeg',sprintf('%s_snr%d.jpg',desc,snr))
% pause
% end%snr   
