function main(opts)

%% opts
% inputs
im_dir = opts.im_dir;
im_fn = opts.im_file_name;

% camera params
focal_len = opts.focal_len;
exposure_time = opts.exposure_time;
readout_time = opts.readout_time;
im_gamma = opts.gamma; % image gamma

% coeff initialization related -- blind deblurring within row segments
nseg = opts.nseg;
iskernelavailable = opts.iskernelavailable;
isusingjiakernel = opts.isusingjiakernel;
blur_kernel_size = opts.blur_kernel_size; 

% coeff initialization related -- poly fitting
psf_bw_thres = opts.psf_bw_thres;
psf_trace_dir = opts.psf_trace_dir;
poly_N = opts.poly_N;
coeffmulfac = opts.coeffmulfac; % slightly increase kernel size, similar to dilatation

% for kernel visualization
isfirstcoeffinit = opts.isfirstcoeffinit; % for getting a glimpse of kernel before proceeding
kgamma = opts.kgamma; 
kgap = opts.kgap;

% for creating Kmat
nsmp = opts.nsmp;

% iteration related params
max_iter = opts.max_iter;
mu = opts.mu;
mu_lastiter = opts.mu_lastiter;
mu_final = opts.mu_final;

% image padding and masks
padding = opts.padding;
padding_final = opts.padding_final;
issaturaredmask = opts.issaturaredmask;

if ~isfirstcoeffinit
    testdir = sprintf('%s/%s',im_dir,opts.testsubdir);
    mkdir(sprintf('%s',testdir));
    save(sprintf('%s/opts.mat',testdir),'opts');
    logFileID = fopen([im_dir '/' opts.testsubdir '/log.txt'], 'w');
end



%% pre-processing
disp(im_dir);

t_exposure = exposure_time;
% t_readout = 1/frame_rate*vertical_resolution_N/vertical_resolution_S * usertune;
t_readout = readout_time;
% normalize so that t_e + t_r = 1
t_exposure = t_exposure/(t_exposure+t_readout);
t_readout = 1-t_exposure;

poly_type = sprintf('poly%d',poly_N-1);
if poly_N == 3
    poly_func = @(coeff,t) coeff(1)*t.^2+coeff(2)*t+coeff(3);
elseif poly_N == 4
    poly_func = @(coeff,t) coeff(1)*t.^3+coeff(2)*t.^2+coeff(3)*t+coeff(4);
elseif poly_N == 5
    poly_func = @(coeff,t) coeff(1)*t.^4+coeff(2)*t.^3+coeff(3)*t.^2+coeff(4)*t+coeff(5);
else
    error('modify poly_func');
end

bc = im2double(imread(sprintf('%s/%s.png',im_dir,im_fn)));
b = rgb2gray(bc);
[h,w,~] = size(b);

K = build_intrisic_mat(size(b),size(b),focal_len);
b = b.^im_gamma;
bc = bc.^im_gamma;



%% fit global camera motion polynomials
% trace each seg
[imseg,~] = gen_im_segs(bc,nseg,im_dir,im_fn);
ker = cell(nseg);
if ~iskernelavailable
    for segindh = 1:nseg(1)
        ker{segindh} = single_level(im_dir,im_fn, imseg{segindh}, 1, segindh, blur_kernel_size);
    end
    save(sprintf('%s/%s_su_kernels.mat',im_dir,im_fn),'ker');
end

% trace
paramk = cell(nseg);
segh = floor(h/nseg(1));
for segindh = 1:nseg(1)
    ts = ((segindh-1)*segh) / (h-1) * t_readout;
    te = min( 1, (segindh*segh-1)/(h-1)*t_readout + t_exposure );
    
    % whether using alternative kernels, in this case, Xu and Jia, [17] in paper
    if isusingjiakernel
        tmp = load(sprintf('%s/%s_jia_kernels.mat',im_dir,im_fn));
    else
        tmp = load(sprintf('%s/%s_su_kernels.mat',im_dir,im_fn));
    end
        
    ker = tmp.ker{segindh};
    paramk{segindh} = psf_traverse(ker,psf_bw_thres(segindh),psf_trace_dir,ts,te,isfirstcoeffinit);
end

% align
t_ref = linspace(0,1,100);
t_patch = cell(nseg(1),1);
y_patch = cell(nseg(1),1);
x_patch = cell(nseg(1),1);
for i = 1:nseg(1)
    t_patch{i} = paramk{i}.end2end(:,4);
    y_patch{i} = paramk{i}.end2end(:,1);
    x_patch{i} = paramk{i}.end2end(:,2);
end
offset_y = zeros(size(nseg));
offset_x = zeros(size(nseg));
for i = 2:nseg(1)
    % locate the overlap where aligement can be based on
    t_overlap_l = t_patch{i}(1);
    t_overlap_r = t_patch{i-1}(end);
    tmp = abs(t_overlap_l - t_ref); [~, ind_l] = min(tmp);
    tmp = abs(t_overlap_r - t_ref); [~, ind_r] = min(tmp);
    t_ref_overlap = t_ref(ind_l:ind_r);
    yoftr = fit(t_patch{i}, y_patch{i}, poly_type);
    y_polyr = poly_func(coeffvalues(yoftr), t_ref_overlap);
    yoftl = fit(t_patch{i-1}, y_patch{i-1}, poly_type);
    y_polyl = poly_func(coeffvalues(yoftl), t_ref_overlap);
    offset_y(i) = mean(y_polyl-y_polyr) + offset_y(i-1);
    xoftr = fit(t_patch{i}, x_patch{i}, poly_type);
    x_polyr = poly_func(coeffvalues(xoftr), t_ref_overlap);
    xoftl = fit(t_patch{i-1}, x_patch{i-1}, poly_type);
    x_polyl = poly_func(coeffvalues(xoftl), t_ref_overlap);
    offset_x(i) = mean(x_polyl-x_polyr) + offset_x(i-1);
end
% align, get the init coeff
figure(1);
t = paramk{1}.end2end(:,4); 
y = paramk{1}.end2end(:,1);
for i = 2:nseg(1)
    t = [t; paramk{i}.end2end(:,4)];
    y = [y; paramk{i}.end2end(:,1) + offset_y(i)];
end
% y = y./K(1);
yoft = fit(t, y, poly_type);
plot(yoft, t, y);
hold on;
% fit ... again for trans_x
t = paramk{1}.end2end(:,4); 
x = paramk{1}.end2end(:,2);
for i = 1:nseg(1)
    t = [t; paramk{i}.end2end(:,4)];
    x = [x; paramk{i}.end2end(:,2) + offset_x(i)];
end
% x = x./K(1);
xoft = fit(t, x, poly_type);
plot(xoft, t, x);
hold off;
saveas(1, sprintf('%s/%s_trans_xy_init.png',im_dir,im_fn), 'png');
close(1);



%% get the coeff from cfit -- this is for the t_x and t_y
coeff = [coeffvalues(yoft)./K(1) -coeffvalues(xoft)./K(1)];
coeff = center_coeff(coeff,poly_func,poly_N);
if isfirstcoeffinit
    Kmat = gen_conv_mat(h,w,1,h, K,coeff,poly_N,poly_func, t_exposure,t_readout, nsmp);
    k = vis_camera_motion_traj_as_psfs(h,w,kgap,Kmat);
    k_output = k;
    imwrite(k_output.^kgamma,sprintf('%s/%s_k_init.png',im_dir,im_fn));
    return
end
coeff = coeff.*coeffmulfac;

% pad to reduce boundary artifacts during the alternating opt part
bp = padarray(b,[padding padding],'replicate');
K = build_intrisic_mat(size(b),size(bp),focal_len);
[hp,wp] = size(bp);
if issaturaredmask
    maskim = zeros(size(bp));
    maskim(bp>250/255) = 1;
    maskim = imfilter(maskim,fspecial('gaussian'));
    maskim(maskim>0) = 1;
    maskim = abs(maskim-1);
else
    maskim = ones(size(bp));
end
mask = spdiags(maskim(:), 0, hp*wp, hp*wp);
% D matrix
Dx = spdiags(ones(hp*wp,1), 0, hp*wp, hp*wp)...
    -spdiags(ones(hp*wp,1), -hp, hp*wp, hp*wp)...
    -spdiags(ones(hp*wp,1), hp*wp-hp, hp*wp, hp*wp);
Dy = spdiags(ones(hp*wp,1), 0, hp*wp, hp*wp)...
    -spdiags(ones(hp*wp,1), 1, hp*wp, hp*wp)...
    -spdiags(ones(hp*wp,1), 1-hp*wp, hp*wp, hp*wp);
D = [Dx; Dy];



%% core rsmb part    
coeff_history = zeros(max_iter,numel(coeff));
mu_mulfac = (mu_lastiter/mu)^(1/(max_iter-1));
for iter = 1:max_iter
    fprintf(logFileID, 'iter: %d',iter);
    fprintf(logFileID, '  %.4f',coeff); fprintf(logFileID, '\n');

    Kmat = gen_conv_mat(hp,wp,1-padding,hp-2*padding,...
        K,coeff,poly_N,poly_func,...
        t_exposure,t_readout, nsmp);
    k = vis_camera_motion_traj_as_psfs(hp,wp,kgap,Kmat);
    k_output = k(padding+1:end-padding,padding+1:end-padding);
    imwrite(k_output.^kgamma,sprintf('%s/%s_k_%d.png',testdir,im_fn,iter));

    if iter == 1
        l = bp;
        l_output = l;
    end
    Klbp = reshape(Kmat*l(:)-bp(:),hp,wp);
    Klb = Klbp(padding+1:end-padding,padding+1:end-padding);
    objval1 = 1/2*norm(Klb(:))^2;
    [lgx, lgy] = imgradientxy(l_output);
    objval2 = mu*(norm(lgx(:),1)+norm(lgy(:),1));
    fprintf(logFileID, 'after K: obj1: %f, obj2: %f, obj: %f, mu: %f\n',...
        objval1, objval2, objval1+objval2, mu);

    l = admmmat(bp,Kmat,D,mu,mask);
    l_output = l(padding+1:end-padding,padding+1:end-padding);
    imwrite(l_output.^(1/im_gamma),sprintf('%s/%s_l_%d.png',testdir,im_fn,iter));

    Klbp = reshape(Kmat*l(:)-bp(:),hp,wp);
    Klb = Klbp(padding+1:end-padding,padding+1:end-padding);
    objval1 = 1/2*norm(Klb(:))^2;
    [lgx, lgy] = imgradientxy(l_output);
    objval2 = mu*(norm(lgx(:),1)+norm(lgy(:),1));
    fprintf(logFileID, 'after l: obj1: %f, obj2: %f, obj: %f, mu: %f\n',...
        objval1, objval2, objval1+objval2, mu);

%     [bx,by] = gradient(bp,1);
%     [lx,ly] = gradient(l,1);

    coeff_history(iter,:) = coeff;

    if iter ~= max_iter
        for iii = 1:1
            coeff = funcGN(coeff,poly_N,poly_func,hp,wp,l(:),bp(:),...
                1-padding,hp-2*padding,K,nsmp,t_exposure,t_readout);
%             coeff = funcGN(coeff,poly_N,poly_func,hp,wp,[lx(:),ly(:)],[bx(:),by(:)]...
%                 ,1-padding,hp-2*padding,K,nsmp,t_exposure,t_readout);
        end
        coeff = center_coeff(coeff,poly_func,poly_N);
    end
    mu = max(mu_lastiter,mu*mu_mulfac);
end
imwrite(l,sprintf('%s/l.png',testdir));
save(sprintf('%s/coeff_history.mat',testdir),'coeff_history');
fclose(logFileID);



%% final
% pad more for the final deblurring
bcp = padarray(bc,[padding_final padding_final],'replicate');
K = build_intrisic_mat(size(b),[size(bcp,1),size(bcp,2)],focal_len);
Kmat = gen_conv_mat(h+2*padding_final,w+2*padding_final,1-padding_final,h,...
    K,coeff,poly_N,poly_func, t_exposure,t_readout, nsmp);
k = vis_camera_motion_traj_as_psfs(h+2*padding_final,w+2*padding_final,kgap,Kmat);
k_output = k(padding_final+1:end-padding_final,padding_final+1:end-padding_final);
imwrite(k_output.^kgamma,sprintf('%s/%s_final_k.png',testdir,im_fn));

if padding_final-padding~=0
    maskim = padarray(maskim,[padding_final-padding padding_final-padding],'replicate');
    mask = spdiags(maskim(:), 0, numel(maskim), numel(maskim));
    % D matrix
    [hpf,wpf,~] = size(bcp);
    Dx = spdiags(ones(hpf*wpf,1), 0, hpf*wpf, hpf*wpf)...
        -spdiags(ones(hpf*wpf,1), -hpf, hpf*wpf, hpf*wpf)...
        -spdiags(ones(hpf*wpf,1), hpf*wpf-hpf, hpf*wpf, hpf*wpf);
    Dy = spdiags(ones(hpf*wpf,1), 0, hpf*wpf, hpf*wpf)...
        -spdiags(ones(hpf*wpf,1), 1, hpf*wpf, hpf*wpf)...
        -spdiags(ones(hpf*wpf,1), 1-hpf*wpf, hpf*wpf, hpf*wpf);
    D = [Dx; Dy];
end

lc = cell(3,1);

for ic = 1:3
    lc{ic} = admmmat(bcp(:,:,ic),Kmat,D,mu_final,mask);
end

l = zeros(h+2*padding_final,w+2*padding_final,3);
l(:,:,1) = lc{1}; l(:,:,2) = lc{2}; l(:,:,3) = lc{3};
l_output = l(padding_final+1:end-padding_final,padding_final+1:end-padding_final,:);
imwrite(l_output.^(1/im_gamma),sprintf('%s/%s_final_l.png',testdir,im_fn));


