%% 时间序列图像的输入
%目前只能输入100以内的帧数，并且最后一位是0，包括10,20,30,40,50,60,70,80,90.
%如果你的时间序列不是这个数目，可以用ImageJ删除一些做伪彩用不到的帧数。
clc; clear all; close all;
%输入多少帧，在下面的变量num,num1中进行定义，比如如果输入70帧，定义下面的num=9,
% 而num1=6
num = 9;  
num1 = 2;  %可修改
Img = [];
      % 输入的图片序列需要用ImageJ 转换成8 bit格式。然后修改下面的路径： 
      % E:\\wuling-weicai\\ct2\\ct00%d%d.tif 作为输入
  for j = 0 : num1
    for i = 0 : num
     filename = sprintf('E:\\image\\20161128 PKA PKC\\293t M3-BFP2 pH\\forskolin 8h\\5 3 cells\\tu\\a00%d%d.tif',j, i);  % 可修改
    I = imread(filename);
    Img(:, :,j*10+i+1) = I;
    end
  end


  
%%  计算F0,以及F0的方差
%   变量num2 定义计算F0需要平均的帧数，根据自己实验进行修改.比如需要平均5帧，则num2=4.
%   变量stdevimage： 表示F0的方差。
num2 = 5;      %可修改                
I = double(I);                 
imsize = size(I);  
sumimage = zeros(imsize); 
sum2image = zeros(imsize);
varianceimage = zeros(imsize);
stddevimage = zeros(imsize);
Fo = zeros(imsize);

% 如果F0不是从第一帧开始，下面的两个i+1,可以自行修改，比如是从第5帧到第10帧，i+1
%修改为i+5,num2等于5.
for i = 0 : num2
    sumimage(:,:) = sumimage(:,:) + double(Img(:,:,i+1));  %可修改
    sum2image(:,:) = sum2image(:,:) + double(Img(:,:,i+1)).^2; %可修改
end
Fo(:,:) = sumimage(:,:)/(num2+1);  % 计算得到平均后的F0,
varianceimage(:,:) = (sum2image(:,:)/(num2+1)) - Fo(:,:).^2;
stddevimage(:,:) = sqrt(varianceimage(:,:));   % 计算得到平均后的F0的方差

 

%%  计算得到信号值： deltaF/F0
deltaFimageseries = [];
F = zeros(imsize);

% 得到一个每一帧减掉F0后的时间序列图像矩阵,即deltaF的矩阵detaFimageseries。并且取阈值。如果像素点的值小于1倍（或2倍，0.5倍等等，根据自己的实验可以修改，最后写文章说明即可，一般都是取了1倍或2倍）方差，
% 将该像素点的值取为0值。
 for j = 0 : num1
    for i = 0 : num
     for n = 1 : imsize(1)
       for m = 1 : imsize(2)
           deltaFimageseries(n, m,j*10+i+1) = Img(n, m,j*10+i+1)-Fo(n,m);
               if ( abs(deltaFimageseries(n,m,j*10+i+1))<=2*stddevimage(n,m))  %可修改,比如改成2*stddevimage(n,m)
                   deltaFimageseries(n,m,j*10+i+1) = 0;
               else
                   deltaFimageseries(n,m,j*10+i+1) = deltaFimageseries(n,m,j*10+i+1);
               end
       end
     end   
    end
 end
 
 
%F0取阈值 ，计算deltaF/F0， 得到一个deltaF/F0的一个时间序列图像矩阵ratioimageseries.
ratioimageseries = [];
for j = 0 : num1
    for i = 0 : num
     for n = 1 : imsize(1)
       for m = 1 : imsize(2)
           if (Fo(n,m)<=0.5)
                   ratioimageseries(n,m,j*10+i+1) = 0;
           else
                   ratioimageseries(n,m,j*10+i+1) = deltaFimageseries(n,m,j*10+i+1)./Fo(n,m);
           end
       end
     end   
    end
end



%%  作伪彩图
  averagefilm(:,:,:)=double(ratioimageseries(:,:,:));
  
%对矩阵进行平滑滤波除噪
for j = 0 : num1
    for i = 0 : num
       averagefilm(:,:,10*j+i+1) = medfilt2(averagefilm(:,:,10*j+i+1 ),[4,4]);
    end
end

%输出伪彩moive
for j = 0 : num1
    for i = 0 : num
figure(1),imshow(averagefilm(:,:,10*j+i+1),[0,0.4],'Border','tight'); % 伪彩的范围[0,3]可修改
colormap(jet);
set(gcf,'visible','on');
str=['E:\image\20161128 PKA PKC\293t M3-BFP2 pH\forskolin 8h\5 3 cells\wc\' num2str(i+10*j),'.bmp'];%输出地路径可以修改
print(1,'-dbmp','-r32',str); 
    end
end





  
