Single_Focus_Confocal_Tracking/Main_code.m at main · LabPresse/Single_Focus_Confocal_Tracking · 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
%%%%%%%%%%%%%%% Single molecule tracking analysis v1.1 %%%%%%%%%%%%%%%%%%%%%
%%%% In the case of problem, please contact
%%%% Sina Jazani (sjazani@asu.ed) or Dr. Steve Presse (spresse@asu.edu).
%%%% The code written by Sina Jazani.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 0
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Clear the history %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
clear
close all

% Add the function folder to the path
addpath('Functions')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 1
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Deltat      =         0.0001   ;     % time interval (Known)
mu          =    100000        ;     % Intensity or average number of photons which came from a particle (Known)
mu_back     =      1000        ;     % Intensity or number of photons in the background (Known)
Diffusion   =         10       ;     % Diffusion coefficent of the molecule
L_xy        =         1        ;     % Periodic boundary in x and y axis
L_z         =         1        ;     % Periodic boundary in z axis
wxy         =         0.3      ;     % Semi-axis of confocal in xy plane (micro meter)
wz          =         0.7      ;     % Semi-axis of confocal in z axis   (micro meter)
Number_step =      1000        ;     % Number of steps

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 2
%%%%%%%%%%%%%%%%%%%%% Generate the synthetic trace %%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[X_real,Y_real,Z_real,Signal] = Sample_Generator( ...
...
...
L_xy  ,L_z, Number_step ,  wxy,wz, Diffusion , Deltat , mu  , mu_back );


% Plot the generated trace
plot( (0:1:Number_step).*Deltat , Signal)
xlim([0 Number_step*Deltat])
xlabel('Time (s)')
ylabel({'Intensity','Number of observed photons per time windows'})

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 5
%%%%%%%%%%%%%%% Define parameters of the model %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Run the Gibbs sampler function based on the number of the iteration.
% There are no automatic MCMC termination for this case and the user
% needs to choose a large number of iteration to target the correct
% posterior distribution.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Please adjust the prior parameters according to the graund true.
%%% If the prior be so off from the actual value, and the likelihood be
%%% weak compare to the prior, the posterior will follow the pior and
%%% sampler will not work correctly
D_alpha_prior  =          5        ;     % Alpha parameter of inverse gamma distr. as the prior of Diffusion coefficient
D_beta_prior   =          50       ;     % Beta parameter of inverse gamma ditsr. as the prior of Diffusion coefficient
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu_prior_xy    =          0.1      ;     % Mean value for X0 prior in x and y coordintes
mu_prior_z     =          0.1      ;     % Mean value for X0 prior in z coordinte
var_prior_xy   =          1        ;     % Variance for X0 prior in x and y coordintes
var_prior_z    =          1        ;     % Variance for X0 prior in z coordinte

%%%%%%%%%%%%%%%%%%%%%%%%%%  Initial vslues  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D              =          10        ;     % Initial Value for Diffusion coefficient
x              = abs(0.01*randn(3,length(Signal))) ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%  Filter type  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Filter_type    =          2        ;     % Filter type (UKF -> Filter_type = 1 and EKF -> Filter_type = 2)
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 6
%%%%%%%%%%%%%%%%%%%%%%%%% Run the algorithm %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Since this process is computationaly heavy, we strongly suggest to
% measure the running time with low number of iterations first.


Max_iter          =      1000;    % Number of iteration

tic
[ D   , x   ] = Gibbs_sampler(...
...
Signal     , Deltat       , D_alpha_prior  , D_beta_prior    , ...
Max_iter   , wxy          , wz             , mu_prior_xy     , ...
mu_prior_z , var_prior_xy , var_prior_z    , mu_back         , ...
x          , D            , mu             , length(Signal)  , ...
Filter_type);
toc
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 7
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Show results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Plot_data( x  , D, X_real,Y_real,Z_real,Signal,wxy,wz,Diffusion)

%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Step 8
%%%%%%%%%%%%%%%%%%%%%%%%%%%% Save results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
save('Results','D','Signal','Deltat','wxy','wz','mu_prior_xy',...
    'mu_prior_z','var_prior_xy','var_prior_z','D_alpha_prior',...
    'D_beta_prior','X_real','Y_real','Z_real','Diffusion','mu','mu_back','x')