function calibrate_chshift(fname,chshift_fname,ch_seq,transformtype)
% take the tetraspeck images in nd2 format, max project in imagej, then save the max
% projection in a multi-tiff. The channels should be in the order of G/T/A/C

if ~exist('ch_seq','var')
    ch_seq=1:4;
end

if ~exist('transformtype','var')
    transformtype='affine';
end


%%
im=imread(fname,1);
for n=2:4
    im(:,:,n)=imread(fname,n);
end

im=im(:,:,ch_seq); %fixing channel sequence


%normalize images
%im=double(im)./max(double(im),[],3);
%% register all channels to G
rng(1);
chshift20x_tform={};
[optimizer,metric]=imregconfig('multimodal');
optimizer.InitialRadius=optimizer.InitialRadius/20;
optimizer.GrowthFactor=1.001;
optimizer.Epsilon=1e-9;
optimizer.MaximumIterations=3000;


parfor n=1:4
   chshift20x_tform{n}=imregtform(im(:,:,n),im(:,:,1),transformtype,optimizer,metric);
end

%%

im_tform=im;

for n=1:4
    im_tform(:,:,n)=imwarp(im(:,:,n),chshift20x_tform{n},'OutputView',imref2d(size(im(:,:,1))));
end

im_lim=prctile(im_tform,99.9,[1,2]);
im1=double(im_tform)./double(im_lim);

for n=1:4
    figure;imshowpair(im1(:,:,1),im1(:,:,n));
    title(['channel',num2str(n)]);
end

save(chshift_fname,'chshift20x_tform');

