
N = 4096;
stackname = 'ribo70S';
filename = [stackname '.st'];
outputfilename = [stackname '_iTCTF.st'];
transformfile = [stackname '.xf'];

pixelsize = 0.3748;

params.voltage  = 300e3;
params.C_s      = 2.7e-3;
params.lambda   = wavlen(params.voltage);
params.fov      = pixelsize*1e-9*N;

fid = fopen(transformfile,'r');
transform = fscanf(fid, '%f %f %f %f %f %f',[6, Inf]);
fclose(fid);

tilt_axis_rotation = atan2(transform(2,:),transform(1,:));
tilt_axis_rotation = tilt_axis_rotation';
tilt_axis_rotation = tilt_axis_rotation + pi; %sometimes FEI header is +- 180 degrees!

[~,def,~,alpha,beta] = getStackInfo(stackname);

beta_r = beta - tilt_axis_rotation;

% load the aligned tilt-series
a = tom_mrcread(filename);
a = a.Value;

% a = int16(rebin(a,[2 2 1])/4);
% N = N/2;

nangles = size(a,3);

edge_blur = 25;
windowing = (0.5-0.5*cos((N/2-abs(xx(N,N)))/edge_blur*pi))*(abs(xx(N,N))>N/2-edge_blur)+(abs(xx(N,N))<=N/2-edge_blur);
windowing = windowing*rot90(windowing);

% allocate output memory
a_iTCTF = double(zeros(size(a)));
tic;
for k = 1:nangles
    [k nangles]
    newparams = params;
    newparams.df_estimated = def(k)*1e-9;

    projection = a(:,:,k);
    projection = dip_image(double(projection));

    % masking the edges
    mean_projection = mean(projection);
    projection = (projection-mean_projection)*windowing+mean_projection;

    % correcting difference between dip_image and imod
    projection = rot90(projection);

    % do the CTF correction
    projection = -2*iCTF_Tilt_2d_projection_adaptive(projection,alpha(k),newparams,beta_r(k));
    
    % store result
    a_iTCTF(:,:,k) = double(rot90(projection,-1));
end
toc;

%write result to disk
tom_mrcwrite(a_iTCTF,'name',outputfilename);
