VIPR_dipole/Main_dipole.m at main · Lew-Lab/VIPR_dipole · GitHub

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
106
107
108
109
110
111
112
113
114
115
116
%%
% written by Boris Ferdman
% this functions performs phase retrieval via a direct numerical gradient
% associated paper: VIPR: Vectorial Implementation of Phase Retrieval for
% fast and accurate microscopic pixel-wise pupil estimation 2020

% License to use this code is granted to all interested, as long as the original author is
% referenced as such. The original author maintains the right to be solely associated with this work.
% Copyright by Boris Ferdman: borisferd@gmail.com
%

clear all;close all;clc
%add current folders file into path
%% data set
VIPR_user_input_dipole;

%%
% don't let plotsize be larger than data
IS.plotsize = min(IS.plotsize,IS.FOV_size);
% define positions per image - [x,y,z,NFP]
q_cord = [z_pos';z_stack_pos']';
IS.psi = q_cord;

%% Gaussian noise estimation per Z
if noisy_flag
    dx=IS.corner_size;
    for j = 1:size(IMG_T,3)
        tmp = [IMG_T(1:dx,1:dx,j);IMG_T(end-dx+1:end-1+1,end-dx+1:end-1+1,j);IMG_T(end-dx+1:end-1+1,1:dx,j);IMG_T(1:dx,end-dx+1:end-1+1,j)];
        %
        mean_stack(j) = mean(tmp(tmp>0));
        % Estimate the std
        std_stack(j) = std(tmp(tmp>0));

    end
    IMG_bu = IMG_T;
else
    IMG_bu = IMG_T;
    std_stack = zeros(1,size(q_cord,1));
    mean_stack = zeros(1,size(q_cord,1));
end
%% thr and reduce offset from the images
IMG_T = IMG_bu;
%
for j = 1:size(IMG_T,3)
    tmp = IMG_T(:,:,j)-mean_stack(j);
    % threshold on %
    if IS.I_thr_flag==1
        maskT = tmp>IS.I_thr.*max(tmp(:));
    else
        maskT = tmp>IS.I_thr.*std_stack(j);
    end
    % erode&dilate mask
    se = strel('disk',1,6);
    erodedI = imerode(double(maskT),se);
    se = strel('disk',1,6);
    mask(:,:,j) = imdilate(erodedI,se);

    IMG_T(:,:,j) = tmp.*mask(:,:,j);

    ff=figure(11)
    subplot(1,3,1)
    imagesc(tmp); colorbar();
    axis  image
    title('input image');
    set(gca,'FontSize',10)
    subplot(1,3,2)
    imagesc(IMG_T(:,:,j)); colorbar();
    axis  image
    title('thresholded image');
    set(gca,'FontSize',10)
    subplot(1,3,3)
    imagesc(IMG_T(:,:,j)-tmp); colorbar();
    axis  image
    title('diff');
    set(gca,'FontSize',10)

end
%close(11)


%% flip the image based on the simdipole code
for j=1:size(IMG_T,3)
    tmp = IMG_T(:,:,j);
    if vec_model_pol=='x'
        tmp = flipud(fliplr(tmp)).';
        %tmp = rot90(rot90(tmp,2),3);

    elseif vec_model_pol == 'y'
        %tmp  = fliplr(flipud(tmp));
        tmp  = flipud(tmp);

    end
    IMG_T(:,:,j)=tmp;
end
%% Fresnel matching of the grid

mask_size = IS.sampling_size;

%% pad data to match simulation size
tmp_stack = [];
for j = 1:size(IMG_T,3)
    tmp = padarray(IMG_T(:,:,j), floor([mask_size - (size(IMG_T,1)), mask_size - (size(IMG_T,2))]/2),'both');
    if mod(mask_size,2)==0
        tmp = padarray(tmp, floor([1,1]),'pre');
    end
    if mod(size(tmp,1),2)==0
        tmp = padarray(tmp, floor([1,1]),'pre');
    end
    % complete 0 to mean noise
    tmp_stack(:,:,j) = tmp;
    tmp = [];
end
% IMG_T = tmp_stack;

%% retrieve the mask
[maskRec,gB,Nph,I_mod] = PR_coverslip_data_v2(IS,tmp_stack,q_cord,std_stack,gpu_flag,vec_model_flag,cost_function_flag,plot_flag,Alg_flag,est_gBlur_flag,noisy_flag,vec_model_pol,initial_pmask_name);