%Analyze a sequence of DIC pics and fluo stacks to quantify beads
%brightness. 
function beads_brightness(path,a)
FR_H=0;
FR_W=0;

manual=0;
if exist([path,'\manual_planes.txt']);
    manual=1;
    plane=load([path,'\manual_planes.txt']);
end
p=0; %number of the bead
scrsz = get(0,'ScreenSize');
bg_sh=10; %defines how far to the left background is sampled
FR=a/2.;
FR1=floor(2*a/6.);%a smaller frame to improve contrast when searching for the plane in focus
k=1;
while k>0
% CHOOSE SINGLE BEADS IN DIC----------------------------------------------     
    fdic=[path,'\','DIC_',num2str(k),'.tif'];
    if exist(fdic)>0
        DIC=imread(fdic);
        f1=figure;
        h1=axes;
        dx=0;
        dy=0;
        imagesc(DIC,[20000 45000]);
        
            line([FR+dx FR+dx],[FR+dy FR_W+dy-FR],'Color','c','LineStyle','-'); %RECTANGLE WITHIN WHICH TO CHOOSE DIC BEADS
            line([FR_H+dx-3*FR-4 FR_H+dx-3*FR-4],[FR+dy FR_W+dy-FR],'Color','c','LineStyle','-');
            line([FR+dx FR_H+dx-3*FR-4],[FR+dy FR+dy],'Color','c','LineStyle','-');
            line([FR+dx FR_H+dx-3*FR-4],[FR_W+dy-FR FR_W+dy-FR],'Color','c','LineStyle','-');
                  
        colormap(gray);
        axis equal
        
       
        but=1;
        j=0;
        while but~=3
            axes(h1);
            [ax,ay,but]=ginput(1);
            if but==3
                break;
            else
                j=j+1;
                axes(h1)
                text(ax,ay,num2str(p+j),'Color','m');
                x(j)=floor(ax-dx);
                y(j)=floor(ay-dy);
            end
        end
        J=j;
        but=1;
        j=0;
        
%CALCULATE BRIGHTNESS Of BEADS AND BACKGROUND FOR STACK IN FOCUS-----------
        GFP=uint32(imread([path,'\','GFP_',num2str(k),'.tif']));
        for j=1:J
            II(j)=0;
            IIbg(j)=0;
            for l=x(j)-FR:x(j)+FR
                for m=y(j)-FR:y(j)+FR
                    II(j)=II(j)+GFP(m,l);
                    IIbg(j)=IIbg(j)+GFP(m,l+2*FR+bg_sh); % GFP function extra argumets define background sampling shift
                end
            end
        end
        BG=mean(IIbg);
        STK=1;
        
% MARK FLUORESCENT REGIONS CORRESPONDING TO SINGLE BEADS
        close(f1);
        f2=figure;
        set(f2,'Position',[1 49 1280 666]);
        h2=axes;
        imagesc(GFP,[20 3500]);
        colormap(gray)
        axis equal
        j=0;
        while j<J
            j=j+1;
            axes(h2)
            hold on;
            p=p+1;
            signal(p)=II(j);
            bcgr(p)=IIbg(j);
            bead(p)=II(j)-IIbg(j);
            stknum(p)=STK;
            xpos(p)=x(j);
            ypos(p)=y(j);
            picnum(p)=k;
            stdbg(p)=std(IIbg);
            line([x(j)-FR x(j)-FR],[y(j)-FR y(j)+FR],'Color','y','LineStyle','-');
            line([x(j)+FR x(j)+FR],[y(j)-FR y(j)+FR],'Color','y','LineStyle','-');
            line([x(j)-FR x(j)+FR],[y(j)-FR y(j)-FR],'Color','y','LineStyle','-');
            line([x(j)-FR x(j)+FR],[y(j)+FR y(j)+FR],'Color','y','LineStyle','-');
            
            line([x(j)+2*FR+bg_sh-FR x(j)+2*FR+bg_sh-FR],[y(j)-FR y(j)+FR],'Color','m','LineStyle','-');
            line([x(j)+2*FR+bg_sh+FR x(j)+2*FR+bg_sh+FR],[y(j)-FR y(j)+FR],'Color','m','LineStyle','-');
            line([x(j)+2*FR+bg_sh-FR x(j)+2*FR+bg_sh+FR],[y(j)-FR y(j)-FR],'Color','m','LineStyle','-');
            line([x(j)+2*FR+bg_sh-FR x(j)+2*FR+bg_sh+FR],[y(j)+FR y(j)+FR],'Color','m','LineStyle','-');
            
        end
        waitforbuttonpress
        close(f2);
    k=k+1;
    else
        break;
    end
    
    % start of data exclusion cycle
    jj=0;
    std_bg_gen=std(bcgr);
    mean_bg_gen=mean(bcgr);
    while jj < p
        jj=jj+1;
        if (bcgr(jj) > (mean_bg_gen + 2*std_bg_gen)) %I still have my doubts as to what the criterion here should be. This one works and is "safe", but it may fail to discard even quite obviously erratic data.
            exclusion_flag(jj)=1;
        else
            exclusion_flag(jj)=0;
        end
    end 
    nn=0;
    mm=0;
    e_p=0;
    while nn < p
        nn=nn+1;
        if exclusion_flag(nn)==0
            mm=mm+1;
            e_signal(mm)=signal(nn);
            e_bcgr(mm)=bcgr(nn);
            e_bead(mm)=bead(nn);
            e_stknum(mm)=stknum(nn);
            e_xpos(mm)=xpos(nn);
            e_ypos(mm)=ypos(nn);
            e_picnum(mm)=picnum(nn);
            e_p=e_p+1;
        end
          
    end    
    % end of data exclusion cycle    

end
       e_AVERAGE_BRIGHTNESS=mean(e_bead)
       e_SEM=std(e_bead)/sqrt(e_p)
       e_N=e_p
       e_pp=[1:e_p];
       e_table=[e_pp' e_picnum' e_stknum' e_xpos' e_ypos' e_signal' e_bcgr' e_bead'];
       dlmwrite([path,'\e_table.xls'],e_table,'newline','pc','delimiter','\t');
  