MATLAB script: first step of calcium imaging data

calcium_ezPlot.m

%% 読み込みファイル (uigetfileを使うと良い)
clear all;
imname = 'Registered_03.tif';

%% Tiff header情報の の読み取り
tObj = Tiff(imname, 'r');
tags = tObj.getTagNames;
info = imfinfo(imname);

tic
frames = zeros(info(1).Width, info(2).Height, length(info));
for i = 1:length(info)
    if mod(i, 100) == 0    
        disp(i)
    end
    frames(:,:,i) = imread(imname, i);
end
toc

% tic
% i = 1;
% while ~tObj.lastDirectory
%     if mod(i, 100) == 0       
%         disp(i)
%     end
%     frames(:,:,i) = tObj.read();
%     tObj.nextDirectory()
%     i = i+1;
% end
% tObj.close;
% toc
tObj.close;


%% 簡易動画表示
figure; 
imagesc(frames(:,:,1));
axis image;
for i = 1:1371
    imagesc(frames(:,:,i));
    title(num2str(i))
    axis image
    pause(1/30) 
end


%% quick plot (まずデータをながめる)
figure;
imagesc(mean(frames, 3))
colormap gray
axis image

rect = imrect;
roiPos = rect.getPosition;
roiPos = int16(roiPos);
close gcf


% rectangle 部分の輝度をもってくる
roiPosX  = [roiPos(1):roiPos(1)+20];
roiPosY  = [roiPos(2):roiPos(2)+20];

roiMask = logical(frames(:,:,1)*0);
roiMask(roiPosX, roiPosY) = 1;


for i = 1:1371
    masked = roiMask.* frames(:,:,i);
    tmp = masked(roiPosX, roiPos);
    df(i) = sum(tmp(:));
end

figure;
hold on;
plot(zscore(df))
% y2 = conv(df, [1 1 1 1 1 1 1])/7;
y2 = zscore(detrend(df));
plot(y2, 'r')


%% ImageJ で作ったROIの読み込み
% roiRoot = '~/git/RoiSet/';
roiRoot = uigetdir();
files = dir([roiRoot filesep '*.roi']);
for c = 1:length(files)
    rois{c} = ReadImageJROI([roiRoot filesep files(c).name]);
end
roiMask = logical(frames(:,:,1)*0);
% ROIは通常揺れ補正済みの動画の平均像を使用して作る．
% ImageJ でも良いが，家ではざっと半自動でつくったあとGIMPやphotoshopで作ってる


%% マスクの作成
figure;
subplot(1,2,2)
imagesc(mean(frames(:,:,:), 3));
colormap gray; hold on;
axis image
for c = 1:length(files)
    roitmp = rois{c}.mnCoordinates;
    roiMask_eachCell{c} = roiMask*0;
    for p = 1:size(rois{c}.mnCoordinates, 1)
        roiMask(roitmp(p,1), roitmp(p,2)) = 1;
        roiMask_eachCell{c}(roitmp(p,1), roitmp(p,2)) = 1;       
    end
    roiMask_eachCell{c} = imfill(roiMask_eachCell{c});       

    plot(roitmp(:,1), roitmp(:,2))
    text(mean(roitmp(:,1))+5, mean(roitmp(:,2))+5, num2str(c), 'Color', [1 1 1])
end
subplot(1,2,1)
imagesc(mean(frames(:,:,:), 3))
axis image;


%% ROI の一応チェック
% figure;
% for c = 1:109
%     imagesc(roiMask_eachCell{c})
%     axis image;
%     title(['ROI - ' num2str(c)])
%     pause(.1)
% end


%%
resp = zeros(109, 1371); 

% ↑本当はちゃんとこういう定数は最初のほうで決めておく
 
% ROI部分の明るさを取ってくるやつ．
for c = 1:109
    disp(c)
    for i = 1:1371        
        masked = roiMask_eachCell{c}.* frames(:,:,i);
        resp(c,i) = sum(masked(:))/sum(roiMask_eachCell{c}(:));
    end
end
% ↑本当はこの for ループを工夫すると10倍くらい早くなるので頑張って


%%
figure; hold on;

% この帯はデモとして適当に引いたものなので，確認して適切な刺激タイミング位置に描画すること
for i = 1:8
    rectangle('Position', [150*(i-1)+75 0 75 30], 'FaceColor', [.9 0.9 0.9], 'EdgeColor', [1 .9 .9])
end

% plot部分 (ループ使わなくても書けるが速度大差ないためこっちをよく使う)
for c = 1:109
    resp_norm(c,:) = resp(c,:)./mean(resp(c,:));
    resp_filt(c,:) = conv(resp_norm(c,:), [  1 1 1  ])/3;
    plot(resp_norm(c,:) + c*0.5, 'Color',[.9 .5 .5] );
    plot(resp_filt(c,2:end-2) + c*.5, 'k');
end