-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmain_MP_SIM.m
106 lines (90 loc) · 3.54 KB
/
main_MP_SIM.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
% 2D SIM reconstruction
% This file is design to process Multi-plane data as generated by a custom
% multi-plane prism (see https://doi.org/10.1101/773440 and
% https://doi.org/10.1038/s41566-018-0109-4)
addpath(['.',filesep,'funcs'])
cdir = cd;
[fname,fpath] = uigetfile({'*.*';'*.All'}); cd(cdir);
% clear workspace
clear simDec; clear wfDec; clear dataDec;
im = loadData([fpath,fname]);
%% reorder raw camera frame into x,y,z,t stack
disp('Data reordering')
data = reorderPlanes(im);
wf = zeros(size(data,1),size(data,2),8);
for k = 1:8
wf(:,:,k) = mean(data(:,:,:,k),3);
end
%% estimate the transformation and prism transmission
disp('Estimate T(z)')
v1x = 50; v1y = 50; v2x = 450; v2y = 450;
[dx,dy,T,cx,x,cy,y] = estimCoregData(wf,0,v1x,v2x,v1y,v2y);
%% coregister and crop the data
disp('Coregister data')
data = coregData(data,dx,dy,T,cx,x,cy,y);
%% estimate the transformation using beads calibration
disp('Beads calibration calculation')
cal = runcalibLabview;
%% coregister and crop the data
disp('Coregister data using beads calibration')
data = coregDataBeads(data,cal);
%% load coregistered data
[im,fname] = loadData;
temp = reshape(im,[size(im,1),size(im,2),size(im,3)/8,8]);
data = zeros(size(temp)); data(:,:,size(data,3)+2,:) = 0;
data(:,:,2:end-1,:) = temp;
ind = strfind(fname,filesep);
fpath = fname(1:ind(end));
fname = fname(ind(end)+1:end);
%% deconvolve data
nIter = 15;
dataDec = deconvolveData(double(data),nIter);
plotStack(dataDec(:,:,7,:),11)
%%
% compute 2D SIM plane per plane
doMask = 1; doApod = 0.05; ki = 0.44;
if exist('dataDec','var')
[wf,sim,wfDec,simDec] = batchSIM(data,doMask,doApod,dataDec,ki);
else
[wf,sim] = batchSIM(data,doMask,doApod);
end
plotStack(sim(:,:,:,1))
if exist('dataDec','var')
plotStack(simDec(:,:,:,1))
end
pps = 108;
ppsSIM = pps/(size(sim,1)/size(im,1));
ind = strfind(fname,'.');
fnameC = fname(1:ind(end)-1);
%% resample SIM images
if exist('simDec','var')
[sim,simDec] = resampleSIM(sim,sim(:,:,5,1),simDec);
else
sim = resampleSIM(sim,sim(:,:,5,1));
end
ppsSIM = pps/(size(sim,1)/size(wf,1))
% ppsSIM = pps/(size(sim,2)/size(wf,2))
%% saving
writeTIFF(reshape(wf,[size(wf,1),size(wf,2),size(wf,3)*size(wf,4)]),[fpath,fnameC,'_wfMean_',getID(5)],[1/pps 1/pps]);
writeTIFF(reshape(sim,[size(sim,1),size(sim,2),size(sim,3)*size(sim,4)]),[fpath,fnameC,'_SIM',getID(5)],[1/ppsSIM 1/ppsSIM]);
if exist('simDec','var')
writeTIFF (reshape(wfDec,[size(wf,1),size(wf,2),size(wf,3)*size(wf,4)]),[fpath,fnameC,'_wfDec_',getID(5)],[1/pps 1/pps]);
writeTIFF(reshape(simDec,[size(sim,1),size(sim,2),size(sim,3)*size(sim,4)]),[fpath,fnameC,'_SIMDec',getID(5)],[1/ppsSIM 1/ppsSIM]);
end
%% compute and save RBGproj
r = [1 1 1 1]; [2 3 2 3];
if exist('dataDec','var')
saveRGBProj(wf,sim,fpath,fnameC,r,wfDec,simDec)
else
saveRGBProj(wf,sim,fpath,fnameC,r)
end
%% compute resolution for a give plane % https://github.com/Ades91/ImDecorr
p = 2;
resWF = 2*pps/getDcorr(gpuArray(apodImRect(wf(1:400,:,p),20)),linspace(0,0.6,50),20)
figure(243639);subplot(221);imagesc(log(abs(fftshift(fftn(fftshift(apodImRect(wf(:,:,p),20)))))+1))
resWFWdec = 2*pps/getDcorr(gpuArray(apodImRect(wfDec(:,:,p),20)))
subplot(222);imagesc(log(abs(fftshift(fftn(fftshift(wfDec(:,:,p)))))+1)); title('wfDec')
resSIM = 2*ppsSIM/getDcorr(gpuArray(apodImRect(sim(:,:,p),20)),linspace(0,1,30),20)
subplot(223);imagesc(log(abs(fftshift(fftn(fftshift(sim(:,:,p)))))+1)); title('sim')
resSIMdec = 2*ppsSIM/getDcorr(gpuArray(apodImRect(simDec(:,:,p),20)))
subplot(224);imagesc(log(abs(fftshift(fftn(fftshift(simDec(:,:,p)))))+1)); title('simDec')