-
Notifications
You must be signed in to change notification settings - Fork 0
/
G2_P_beads.m
105 lines (105 loc) · 4.86 KB
/
G2_P_beads.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
clear all
startup
% it is the emission wavelength
%% also for registration
Img_C = read3dtiff('C:\Users\ca38cin\Documents\beads\res\C1-1742.tif');
ImgCor603 = correctZdrift(Img_C(:,:,0:end),'MyShifts.txt',1,Img_C(:,:,4));
ImgP603 = mean(ImgCor603,[],3);
Img_C = read3dtiff('C:\Users\ca38cin\Documents\beads\res\C2-1742.tif');% max cc (254,256);center(256,256)
ImgCor719 = correctZdrift(Img_C(:,:,0:end),'MyShifts.txt',1,Img_C(:,:,4));
ImgP719 = mean(ImgCor719,[],3);
Im=cat(3,ImgP603,ImgP719);
ImgTes = correctZdrift(Im(:,:,0:end),'MyShifts.txt',1,Im(:,:,0));
ImgP719 = squeeze(ImgTes(:,:,1));
ImgP603 = squeeze(ImgTes(:,:,0));
%% subtract back ground
CH = 639;%561
BK = ImgP719(20:20+64,67:67+64);
BK = dip_image(repmat(double(BK),2));
Img719 = ImgP719(350:350+129,90:90+129)-BK;
if (CH == 561)
BK = ImgP603(20:20+64,67:67+64);
BK = dip_image(repmat(double(BK),2));
Img603 = ImgP603(350:350+129,90:90+129)-BK;
% k=mean(Img719(65:65+24,90:90+24))/mean(Img603(65:65+24,90:90+24));
k=0.216022765385305;
Img = Img603-k*Img719;
else
Img =Img719;
clear Img719
Img=Img-min(Img);
end
%% PSF
PSF1 = 'reuse'
switch PSF1
case 'reuse'
if (CH == 561)% there are 3 pixels missmatching in x direction
h=load('E:\20220105\restoration\G2_P\resample\measured_PSF\TV\10x\561\PSF_561.mat');
h=h.h;
elseif(CH == 639)
h=load('E:\20220105\restoration\G2_P\resample\measured_PSF\TV\10x\639\PSF_639.mat');
h=h.h;
end
case 'exper'
disableCuda()
[myPSF,FWHM,myResiduum,beadsAt,AllParam,subPixelShifts]=ExtractMultiPSF(ImgP,[25,25],[46.0356,46.0356],1,[])
h = extract (myPSF, size(Img));
case 'calcu'
h = kSimPSF( {'lambdaEm',719;'Pi4Em',0;'na',1.2;'ri',1.33;'sX',size(Img,1);'sY',size(Img,2);'sZ',1;'scaleX',46.0356;'scaleY',46.0356;'scaleZ',1;'lambdaEx',CH;'pinhole',1;'confocal',0;'nonorm',0;'Pi4Ex',0;'computeASF',0;'circPol',0;'scalarTheory',0;'o',''});
end
%% only if using an interpolated image as in put
if(0)
mask=rr(size(Img),'freq')<0.5;
disableCuda()
Img =real(ift(extract(ft(DampEdge(Img)).*mask,size(Img).*10)));
h =real(ift(extract(ft(DampEdge(h)).*mask,size(Img))));
Img=Img-min(Img);
end
%% Load EM
Img_EM = read3dtiff('C:\Users\ca38cin\Documents\beads\res\FFT_TV_w2.tif');
edg=1-(Img_EM==0);
% Img_EM = smooth(Img_EM,1);
image_out = watershed(Img_EM,1,10,0);
Img_EM = smooth(Img_EM,1);
Img_EM_up = Img_EM(:,0:495);
[out_up,thres] = threshold(Img_EM_up,'fixed',160);%ROI2
Img_EM_down = Img_EM(:,496:end);
[out_down,thres] = threshold(Img_EM_down,'fixed',170);%ROI2
out= dip_image(zeros(size(Img_EM)));
out(:,0:495)=out_up;
out(:,496:end)= out_down;
Img_EM =(1-out)*(1-image_out).*edg;
%% deconvolution
Regularization='TG'
NIter=300;
lam=7.5e-13;
ep=1e-5;
useCuda=1;
FP='ForcePos'
% FP='ForcePiecewisePos'
% FP='ForceHyperPos'
switch (Regularization)
case 'GG'
[RefImgX,RefImgY]=RefImg_1(Img_EM,0,2);
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{'CLE_GS',{lam RefImgX RefImgY ep};FP,[];'Resample',10;'NormFac',1},[1,1,1],[10 10],[],useCuda)
% res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{'CLE_GS',{lam,RefImgX,RefImgY ep};FP,[];'NormFac',1},[1,1,1],[100 100],[],useCuda)
case 'IG'
res=GenericDeconvolution(LM,h,NIter,'Poisson',[],{'IG',{lam,Img_EM,ep};FP,[]},[1,1,1],[0 0],[],useCuda);
case 'IT'
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{'IG',{lam,Img_EM,ep};'TV',[1e-7,1e-10];FP,[];'Resample',10;'NormFac',1},[1,1,1],[10 10],[],useCuda);
case 'TV'
% res603=GenericDeconvolution(ImgP,h,NIter,'Poisson',[],{Regularization,[lam,ep];'ForcePos',[];'NormFac',1},[1,1,1],[100 100],[],useCuda);
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{Regularization,[lam,ep];FP,[];'Resample',10;'NormFac',1},[1,1,1],[10 10],[],useCuda);
case 'GR'
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{Regularization,lam;FP,[];},[1,1,1],[10 10],[],useCuda);
% res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{Regularization,lam;FP,[];'Resample',10},[1,1,1],[10 10],[],useCuda);
case 'EG'
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{Regularization,{lam,Img_EM,ep};'Resample',10;FP,[];'NormFac',1},[1,1,1],[10 10],[],useCuda);
case 'TK'
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{'CO',lam;'Resample',10;FP,[];'NormFac',1},[1,1,1],[10 10],[],useCuda);
case 'TE'
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{'EG',{1e-8,Img_EM,1e-9};'CO',lam;'Resample',10;FP,[];'NormFac',1},[1,1,1],[10 10],[],useCuda);
case 'TG'
[RefImgX,RefImgY]=RefImg_1(Img_EM,0,2);
res=GenericDeconvolution(Img,h,NIter,'Poisson',[],{'CLE_GS',{1e-15,RefImgX,RefImgY,1e-5};'CO',lam;'Resample',10;FP,[];'NormFac',1},[1,1,1],[10 10],[],useCuda);
end