%FRET_analysis_PC (4/13/10 Kopito Ronen, Nicolas Pirez)
%function out=FRETanalysis4(method,varargin);
%To use this code, there must be no spaces in file or folder names
%Files to analyze must be .tif files 
clear;clf;
method = 'picksize'; 
nargin = 1;

addpath(pwd);

[fname,path]=uigetfile('*.tif','Select TIFF file to open');
currdir=pwd;
eval(['cd ''' path ''';']);

%generate headings for excel file

d=path;
cd(path);
Destination_Folder = d;
d = uint8(d); d = (d==92); d = find(d,2,'last');
folder_name = Destination_Folder(d(1)+1:end-1);
file_name =  [fname(1:end-4)  '_FRET_Data.xls'];
file_name = [Destination_Folder '\' file_name];
ind =1;
Head1={'Background1' ,'YFP' ,'Background2' ,'CFP' ,'FRET'};
xlswrite(file_name, Head1,'range',num2str(ind));


if nargin>1
    strt = round(double(varargin{1}(1)));
    fin = round(double(varargin{1}(2)));
    stack=tiffread2(fname,strt,fin);
else
    stack=tiffread2(fname);
end
fname_mod=fname(1:end-4);

eval(['cd ''' currdir ''';'])

for i=1:length(stack)
    m(i,:,:)=stack(i).data;
end
m=double(m);

% select regions 1 and 2 (waiting on imaging toolbox)
switch method
    case 'mean'
        im1=squeeze(max(m,[],1));
    case 'pick'
        im1=squeeze(max(m,[],1));
    case 'picksize'
        im1=squeeze(max(m,[],1));
end
figure;
imagesc(im1);
zoom on;
disp('Zoom in ... Press a key to continue to and select REGION in the YFP channel (right)'); pause; zoom off;
[x,y]=getline(gca,'closed');
region1=double(roipoly(im1,x,y));
imagesc(im1);
zoom on;
disp('Zoom in ... Press a key to continue to and select BACKGROUND in the YFP channel (right)'); pause; zoom off;
[x,y]=getline(gca,'closed');
background1=double(roipoly(im1,x,y));
imagesc(im1);
zoom on;
disp('Zoom in ... Press a key to continue to and select REGION in the CFP channel (left)'); pause; zoom off;
[x,y]=getline(gca,'closed');
region2=double(roipoly(im1,x,y));
imagesc(im1);
zoom on;
disp('Zoom in ... Press a key to continue to and select BACKGROUND in the CFP channel (left)'); pause; zoom off;
[x,y]=getline(gca,'closed');
background2=double(roipoly(im1,x,y));   % background subtract!

% integrate over regions
totsize=prod(size(region1));
region1=reshape(region1,1,totsize);
region2=reshape(region2,1,totsize);
background1=reshape(background1,1,totsize);
background2=reshape(background2,1,totsize);

r1=find(double(region1));
r2=find(double(region2));
b1=find(double(background1));
b2=find(double(background2));
for i=1:size(m,1);
    currm=reshape(m(i,:,:),1,totsize);
    switch method
        case 'mean'
            int1(i)=mean(currm(r1));
            int2(i)=mean(currm(r2));
            back1(i)=mean(currm(b1));
            back2(i)=mean(currm(b2));
        case 'pick'
            if i==1
                back1(i)=mean(currm(b1));
                back2(i)=mean(currm(b2));
                hist(currm(r1)-back1(i),50); zoom on;
                p1=input('YFP threshold: ');
                hist(currm(r2)-back2(i),50); zoom on;
                p2=input('CFP threshold: ');
                int1(i)=mean(currm(r1(find(currm(r1)>=back1(i)+p1))));
                int2(i)=mean(currm(r2(find(currm(r2)>=back2(i)+p2))));
            else
                back1(i)=mean(currm(b1));
                back2(i)=mean(currm(b2));
                int1(i)=mean(currm(r1(find(currm(r1)>=back1(i)+p1))));
                int2(i)=mean(currm(r2(find(currm(r2)>=back2(i)+p2))));
            end
        case 'picksize'  % choose brightest pixels
            if i==1
                back1(i)=mean(currm(b1));
                back2(i)=mean(currm(b2));
                imagesc(squeeze(m(1,:,:))); zoom on;
                disp('Zoom in ... Press a key to continue and select ROI in single frame on the YFP channel'); pause; zoom off;
                [x,y]=getline(gca,'closed');
                neursize1=double(roipoly(im1,x,y));
                nsize1=length(find(neursize1));  % just want the number, ignore position
                imagesc(squeeze(m(1,:,:))); zoom on;
%                 disp('choose the single frame neuron outline on the LHS'); pause; zoom off;
%                 [x,y]=getline(gca,'closed');
%                 neursize2=double(roipoly(im1,x,y));
%                 nsize2=length(find(neursize2));  % just want the number, ignore position
                nsize2=nsize1;
                back1(i)=mean(currm(b1));
                back2(i)=mean(currm(b2));
                s1=-1*sort(-1*currm(r1));  % just sort 'em
                s2=-1*sort(-1*currm(r2));
                int1(i)=mean(s1(1:nsize1));
                int2(i)=mean(s2(1:nsize2));
%                 keyboard;
            else
                back1(i)=mean(currm(b1));
                back2(i)=mean(currm(b2));
                s1=-1*sort(-1*currm(r1));  % just sort 'em
                s2=-1*sort(-1*currm(r2));
                int1(i)=mean(s1(1:nsize1));
                int2(i)=mean(s2(1:nsize2));
            end
        otherwise
            p=str2num(method);
            if p<1 & p>0   % percentile not great: as thing gets larger and smaller
                int1(i)=mean(currm(r1(find(currm(r1)>=prctile(currm(r1),p)))));
                int2(i)=mean(currm(r2(find(currm(r2)>=prctile(currm(r2),p)))));
                back1(i)=mean(currm(b1(find(currm(b1)>=prctile(currm(b1),p)))));
                back2(i)=mean(currm(b2(find(currm(b2)>=prctile(currm(b2),p)))));
            elseif p>1 & p>0  % might want one number for each?  how optimize this?
                back1(i)=mean(currm(b1));
                back2(i)=mean(currm(b2));
                int1(i)=mean(currm(r1(find(currm(r1)>=back1(i)+p))));
                int2(i)=mean(currm(r2(find(currm(r2)>=back2(i)+p))));
            elseif p<0
                back1(i)=mean(currm(b1));
                back2(i)=mean(currm(b2));
                int1(i)=mean(currm(r1(find(currm(r1)>=back1(i)*(1-p)))));
                int2(i)=mean(currm(r2(find(currm(r2)>=back2(i)*(1-p)))));
            end
    end            
end
framenum=[1:size(m,1)];
int1=int1-back1;
int2=int2-back2;  % subtract it off...
FRET = int1./int2; %new defined variable YFP/CFP ratio

% plot time series of R1,R2,R1/R2
figure;
subplot(3,2,1);
plot(framenum,back1,'r'); hold on;
plot(framenum,back2,'b');
ylabel('background');
zoom on
subplot(3,2,3);
plot(framenum,int1/max(int1),'r'); hold on;
plot(framenum,int2/max(int2),'b');
ylabel('CFP and YFP');
zoom on
subplot(3,2,5);
plot(framenum,FRET,'k');
ylabel('YFP/CFP');
xlabel('timestep');
set(gca,'ylim',[min(FRET) max(FRET)]); %using new variable FRET
zoom on
subplot(1,2,2);
plot((int1), (int2),'k.-');
xlabel('YFP'); ylabel('CFP');
zoom on

A = int1./int2;
C=corrcoef(diff(int1),diff(int2));
disp(['Correlation between dY and dC is r = ' num2str(C(1,2))]);

% write out to text file with same root name...
Data1=[back1;int1;back2;int2;int1./int2]';
xl_length=length(Data1)+1;
a=['A' num2str(2)];
b=['E' num2str(xl_length)];
xlswrite(file_name, Data1,'range',a);

%writes a dat file in which this variables are saved
file_name = [fname(1:end-4) '_FRET_Analysis.dat'];
file_name = [Destination_Folder '\' file_name];
save (file_name, 'back1' ,'int1' ,'back2' ,'int2' ,'FRET');

%saving the figure
%figure_name = [fname(1:end-4) '_FRET.fig'];
%figure_name = [Destination_Folder '\' figure_name];
%saveas (figure_name, 'figure_name.jpg');

