clearvars
clc
%change microscope name for new scopes!
microscope_name = 'birch';
today = date;

%filename of channel shift matrix
chshift20x_filename = 'chshift20x-birch-20240208.mat';

load(chshift20x_filename)

% which folder contains the calibration files
im_folders={'G_chamber','T_chamber','A_chamber','C_chamber'};

% the sequence of maxproject stack image
channel_seq=[1 2 3 4];%1=G, 2=T, 3=A, 4=C

hybthresh1=[100 50 80 50];

%% shift images

for i=1:numel(im_folders)
    cd(im_folders{i});
    imfiles=dir('MAX_img_stack.tif');
    imfiles={imfiles.name};
    for n=1:numel(imfiles)
        im=imread(imfiles{n});
        for m=2:4
            im(:,:,m)=imread(imfiles{n},m);
            
        end
        ref=imref2d(size(im(:,:,1)));
        t_im=im;
        for m=1:4
            t_im(:,:,m)=imwarp(im(:,:,m),chshift20x_tform{m},'OutputView',ref);
        end
        imwrite(t_im(:,:,1),['shifted',imfiles{n}]);
        for m=2:4
            imwrite(t_im(:,:,m),['shifted',imfiles{n}],'WriteMode','append');
        end

    end
    cd ..
end
%%

%% background subtraction
r=7;%background subtraction radius
se=strel('ball',r,r);
t_im={};
counter=1;
for i=1:numel(im_folders)
    cd(im_folders{i});
    imfiles=dir('shifted*.tif');
    imfiles={imfiles.name};
    for n=1:numel(imfiles)
        im=imread(imfiles{n});
        for m=2:4
            im(:,:,m)=imread(imfiles{n},m);
        end
        t_im{counter}=im;
        
        for m=1:4
            t_im{counter}(:,:,m)=imtophat(im(:,:,m),se);
        end

        imwrite(t_im{counter}(:,:,1),['bgn',imfiles{n}]);
        for m=2:4
            imwrite(t_im{counter}(:,:,m),['bgn',imfiles{n}],'WriteMode','append');
        end
        counter=counter+1;

    end
    cd ..
end

%%
clc

chprofile20x=zeros(4,4);
channel_seq=[1 2 3 4];
for i=1:numel(channel_seq)
    disp(i)
    im=t_im{i};
    r=find_rol(im,hybthresh1);
    sig=zeros(numel(r{channel_seq(i)}),4);
    %find rolonies
    for n=1:4
        im1=im(:,:,n);
        sig(:,n)=im1(r{channel_seq(i)});
    end
    % calculate signal vector for rolonies
    sig_coeff=median(sig(1:100,channel_seq(i)));
    sig=sig./sig(:,channel_seq(i));
    chprofile20x(channel_seq(i),:)=mean(sig);
    %remove noise
    chprofile20x(chprofile20x<0.1)=0;
    %scale by median signal intensity
    chprofile20x(channel_seq(i),:)=chprofile20x(channel_seq(i),:)*sig_coeff;
end

%normalize so that weakest channel has a diagonal of 1.
a=chprofile20x.*eye(4);

chprofile20x=chprofile20x/min(a(a>0));

save(['chprofile20x-', microscope_name,'-',date,'.mat'],'chprofile20x');




%%


function r=find_rol(lim,hybthresh1)
   r={};rsub=r;
    parfor n=1:4
        a=lim(:,:,n);
        CC = bwconncomp(imregionalmax(imreconstruct(max(a-hybthresh1(n),0),a)));
        r{n}=zeros(length(CC.PixelIdxList),1);
        for i=1:length(CC.PixelIdxList)
            [~,I]=max(a(CC.PixelIdxList{i}));
            r{n}(i)=CC.PixelIdxList{i}(I); %linear indexed peak positions
        end
        [y,x]=ind2sub(size(lim(:,:,n)),r{n});
        rsub{n}=[x,y];%rsub is in x, y, consistent with lroi1
    end
end