matlab_foci_GUI/Intensity_weighting.m at master · junlabucsd/matlab_foci_GUI · 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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
clear all; clc;
% close all;

parameters

%% load data
cell_data = load(cell_foci_file_path);

% extract data
L = length(fieldnames(cell_data));
fnames = fieldnames(cell_data);

channel = cell(1,1);
Foci_all = [];


for i = 1:L

    N = length( cell_data.(fnames{i}).times );

    channel = cell_data.(fnames{i});

    if isempty(channel.foci_h) == 0 && channel.division_time

        for j = 1:N

            if isempty(channel.foci_h) == 0

                iscell_foci = double(iscell(channel.foci_h));

                if iscell_foci == 1
                    Foci_all = [Foci_all; [channel.foci_h{1,j}']];
                end

                if iscell_foci == 0
                    Foci_all = [Foci_all; [channel.foci_h(j,1)']];
                end

            end

        end

    end

    if mod(i,100)==0
        i
    end

end

% save('../../analysis/IW_foci_20171026.mat');


%% plot
% close all;
% load('../../analysis/IW_foci_20171026.mat');

foci_counter = length(Foci_all(:,1));
bin_wid_scale = 500;

colors = [46 49 146;
          28 117 188;
          0 167 157;
          57 181 74;
          141 198 63;
          251 176 64;
          241 90 41;
          239 65 54]/255; %illustrator

positions = [400, 400, 420, 400];


% %-------------Fitted Peak intensity of focus------------
%-------------Total intensity of focus------------
bin_wid4 = 2*mean(Foci_all(:,1))/bin_wid_scale;

fig = figure;
set(fig,'Position',positions(1,:));
hold on;

[f, x] = hist(Foci_all(:,1), (max(Foci_all(:,1))-min(Foci_all(:,1)))/bin_wid4);
h1 = plot(x,f/length(foci_counter));
h1.Color = colors(1,:);
set(h1,'LineWidth',1,'LineStyle','-');


set(gca,'YScale','linear','YTick',[],'YTickLabel',{});

% xlabel('Total intensity of focus (AU)','fontsize',20)
% set(gca,'XScale','linear','XTick',[0 30 60 90 120 150 180 210 240]*1e2,'XTickLabel',{'0','','','','1.2\times10^4','','','','2.4\times10^4'})
set(gca,'XScale','linear','XTick',[0 30 60 90 120 150 180 210 240]*2e2,'XTickLabel',{'','','','','','','','',''})
xlim([0 mean(Foci_all(:,1))+5*std(Foci_all(:,1))])

set(gca,'TickLength',[0.025 0.05],'fontsize',20,'TickDir','out','PlotBoxAspectRatio',[1 1 1])

%% fit the intensity weighting distribution
% close all;
%%------multi-gaussian fitting
[fit1 gof] = fit(x.', f.','gauss2');
% [fit1 gof] = fit(x.', f.','gauss2','Lower',[80 5e3 0 20 1.0e4 0],'Upper',[160 8e3 1e4 80 2.5e4 1e4])

g1 = fit1.a1*exp(-((x-fit1.b1)./fit1.c1).^2);
g2 = fit1.a2*exp(-((x-fit1.b2)./fit1.c2).^2);
% g3 = fit1.a3*exp(-((x-fit1.b3)./fit1.c3).^2);
figure;
hold on;

plot(fit1,x,f);
plot(x,g1,'r--');
plot(x,g2,'b--');
% plot(x,g3,'m--');

pg1 = normcdf(x,fit1.b1,fit1.c1/(2^0.5));
pg2 = normcdf(x,fit1.b2,fit1.c2/(2^0.5));

peak_ind_x = zeros(1,length(x));
thr_x = 0;

for i=1:length(x)
    if min(pg1(i),1-pg1(i)) >= min(pg2(i),1-pg2(i))
        peak_ind_x(i) = 1;
    else
        peak_ind_x(i) = 2;
    end
end

for i=1:length(x)-1
    if peak_ind_x(i+1)-peak_ind_x(i)>0
        thr_x = (x(i+1)+x(i))/2;
    end
end