%This program can  the speed of the calcium propagation
%Author:    Xinwei Gao, Qinchun Guo*
%           Imaging Facility of Chinese Institute for Brain
%           Research,Beijing(CIBR)
%           gaoxinwei@cibr.ac.cn
%           guoqingchun@cibr.ac.cn
%%
%????????.tif
clear all 
clc
roi=4;% the sample rate
%%
[fileName,filePath]=uigetfile('.tif');
info=imfinfo([filePath fileName]);
mImage=info.Height;
nImage=info.Width;
numberImage=length(info);
bitDepth=info.BitDepth;
if bitDepth==24
    Image=zeros(mImage,nImage,3,numberImage);    
    h=waitbar(0,'Loading...');
    for i=1:numberImage
        Image(:,:,:,i)=imread([filePath fileName],i);
        if ~rem(i,20)
            waitbar(i/numberImage);
        end
    end
    Image=reshape(Image(:,:,2,:),mImage,nImage,numberImage);
    close(h); 
else
    Image=zeros(mImage,nImage,numberImage);
    h=waitbar(0,'Loading...');
    for i=1:numberImage
        Image(:,:,i)=imread([filePath fileName],i);
        if ~rem(i,20)
            waitbar(i/numberImage);
        end
    end
    close(h); 
end
% %%
% Image1D=sum(Image,2);
% imagesc(squeeze(Image1D)');
% for i=1:numberImage
%     plot(squeeze(Image1D(:,1,i)));
%     pause(0.2);
% end
% %%
% Imagediff=diff(Image,1,3);
% Imagediff(Imagediff<0)=0;
% for i=1:80
%     imagesc(Imagediff(:,:,i));
%     title(int2str(i));
%     pause(0.2);
% end


%%
F0=mean(Image(:,:,1:10),3);
deltaFF=100*(Image-F0)./F0;
for i=1:numberImage
%     imagesc(deltaFF(:,:,i),[min(min(min(deltaFF))),max(max(max(deltaFF)))]);
%     title(int2str(i));
%     pause(0.2);
    imwrite(uint16(deltaFF(:,:,i)),'myMultipageFile.tif','WriteMode','append')  %myMultipageFilej is the DELTAF/F data
end


%%


% for i=1:20
%    imagesc(ImageDeltF(:,:,i));
%    title(int2str(i));
%     pause(0.5);
% end
% imagesc(deltaFF(:,:,30));
figure(1)
imagesc(Image(:,:,100))
[Xcenter,Ycenter]=getpts(figure(1));  %????????


%%
%draw ROI in the figure

for i=1:roi
    
I=Image(:,:,20); % ??????
figure(2)
imagesc(I);
BW_gfproi=imfreehand(gca);
w=waitforbuttonpress;
if w==0
    disp('Button click')
else
    disp('Key press')
end
BW(:,:,i)=createMask(BW_gfproi);% ??????
% [c,r]=ginput(3); % 20???????20?????????????????
% BW(:,:,i)= roipoly(I,c,r); % ????????
end
for i=1:6
   imagesc(BW(:,:,i));
   title(int2str(i));
    pause(0.5);
end

%%
%GXW version-mean intensity
for j=1:size(BW,3)
for i=1:size(Image,3)
     meansignal(i,j)=(sum(sum(Image(:,:,i).*BW(:,:,j))))/sum(sum(BW(:,:,j)));
    end
end
meansignal=meansignal';
for i=1:roi
    plot(meansignal(i,:));
    title(int2str(i));
    hold on
    pause(0.5);
end
%%
timepoint=zeros(1,roi);
distance=zeros(1,roi);
for i=1:roi
%     i=4
%     i=2;
xp=1:1:size(Image,3);
yp=meansignal(i,:);
% dy=gradient(yp)./gradient(xp);
% dx=1:1:size(Image,3);
dx1=2:1:size(Image,3);
% figure
% plot(dx,dy,'r')
dy1=diff(yp)./diff(xp);

if max(dy1)<200
    timepoint(1,i)=nan;
    distance(1,i)=nan;
else
 
[~,pos]=max(dy1);
timepoint(1,i)=pos;
S=regionprops(BW(:,:,i),'Centroid');
centroids(i,:) = cat(1, S.Centroid);
distance(1,i)=sqrt((centroids(i,1)-Xcenter)^2+(centroids(i,2)-Ycenter)^2);
end
figure(i)
plot(dx1,dy1,'b')
legend('??','??')
hold on 
plot(xp,yp)
end
distance=distance(~isnan(distance));
timepoint=timepoint(~isnan(timepoint));
%%
% figure(5)
% clims = [0 3000];
% subplot(3,3,1)
% imagesc(ImageDeltF(:,:,10),clims)
% title('Frame 10')
% subplot(3,3,2)
% imagesc(ImageDeltF(:,:,15),clims)
% title('Frame 15')
% subplot(3,3,3)
% imagesc(ImageDeltF(:,:,20),clims)
% title('Frame 20')
% subplot(3,3,4)
% imagesc(ImageDeltF(:,:,25),clims)
% title('Frame 25')
% subplot(3,3,5)
% imagesc(ImageDeltF(:,:,30),clims)
% title('Frame 30')
% subplot(3,3,6)
% imagesc(ImageDeltF(:,:,35),clims)
% title('Frame 35')
% subplot(3,3,7)
% imagesc(ImageDeltF(:,:,40),clims)
% title('Frame 40')
% subplot(3,3,8)
% imagesc(ImageDeltF(:,:,45),clims)
% title('Frame 45')
% subplot(3,3,9)
% imagesc(ImageDeltF(:,:,50),clims)
% title('Frame 50')
 cftool(timepoint,distance)


%        
  




















