Skip to content

Commit 08b4868

Browse files
anneliesdewaelanneliesdewael
authored
Merge pull request #15 from quantitativeTEM/with_hmm
Addition of hidden Markov Model analysis to count atoms in time series
2 parents ed79cc2 + 0595a4f commit 08b4868

26 files changed

Lines changed: 859 additions & 16 deletions

‎Examples/example_HMM_input.mat‎

1.75 MB
Binary file not shown.

‎GUI/LeftPanels/Analysis/panelAtomCountHMM.m‎

Lines changed: 32 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,32 @@
1+
function panel = panelAtomCountHMM()
2+
% panelAtomCountHMM - Create panel with buttons for counting atoms using
3+
% hidden Markov models
4+
%
5+
% syntax: panel = panelAtomCountHMM()
6+
% panel - structure containing the necessary information to create a
7+
% panel with buttons in StatSTEM
8+
%
9+
10+
%--------------------------------------------------------------------------
11+
% This file is part of StatSTEM
12+
%
13+
% Copyright: 2019, EMAT, University of Antwerp
14+
% Author: K. H. W. van den Bos, A. De wael
15+
% License: Open Source under GPLv3
16+
% Contact: sandra.vanaert@uantwerpen.be
17+
%--------------------------------------------------------------------------
18+
19+
panel = struct;
20+
21+
% Start by defining the name
22+
panel.name = 'Atom counting - Time series';
23+
24+
panel.row = struct;
25+
%% Define buttons, etc row per row (don't put too much buttons per row for size limitations)
26+
% Define buttons, etc for the first row per option. Each row has a maximum
27+
% width of 162 pixels
28+
29+
panel.row(1).option1 = pushbutton('name','Run HMM','width',162,'func','analyseFactorialHMM','input','inputHMM','output','outputHMM','fitting',1,...
30+
'figEnd','Observations Time Series','figOptEnd',{'Atom Counts Time Series'},'reshowFigEnd',true);
31+
32+

‎GUI/RightPanels/possibleImagesStatSTEM.m‎

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -34,6 +34,11 @@
3434
out.option1 = addFigOpt(out.option1,[char(949),'_yy'],'showStrainEpsYY','strainmapping.eps_yy');
3535
out.option1 = addFigOpt(out.option1,[char(969),'_xy'],'showStrainOmgXY','strainmapping.omg_xy');
3636

37+
% Option1b: Observation Time Series
38+
out.option1b = listImagesStatSTEM('Observations Time Series','showObservation_T','inputHMM.obs_T');
39+
out.option1b = addFigOpt(out.option1b,'Coordinates Time Series','plotCoordinates_T','outputHMM.coordinates_T',true); % Show this option by default
40+
out.option1b = addFigOpt(out.option1b,'Atom Counts Time Series','plotAtomCounts_T','outputHMM.H_viterbi');
41+
3742
% Option2: Model
3843
out.option2 = listImagesStatSTEM('Model','showModel','output.model','input');
3944
out.option2 = addFigOpt(out.option2,'Input coordinates','plotCoordinates','input.coordinates');
@@ -52,6 +57,11 @@
5257
out.option2 = addFigOpt(out.option2,[char(949),'_yy'],'showStrainEpsYY','strainmapping.eps_yy');
5358
out.option2 = addFigOpt(out.option2,[char(969),'_xy'],'showStrainOmgXY','strainmapping.omg_xy');
5459

60+
% Option2b: Model Time Series
61+
out.option2b = listImagesStatSTEM('Models Time Series','showModel_T','inputHMM.model_T');
62+
out.option2b = addFigOpt(out.option2b,'Coordinates Time Series','plotCoordinates_T','outputHMM.coordinates_T',true); % Show this option by default
63+
out.option2b = addFigOpt(out.option2b,'Atom Counts Time Series','plotAtomCounts_T','outputHMM.H_viterbi');
64+
5565
% Option3: Histogram SCS
5666
out.option3 = listImagesStatSTEM('Histogram SCS','showHistogramSCS','output.volumes');
5767
out.option3 = addFigOpt(out.option3,'GMM components','plotGMMcomp','atomcounting.selMin',true);
@@ -89,4 +99,4 @@
8999
out.option12 = addFigOpt(out.option12,'Coordination number','showModelCoorNum','model3D.coorNum');
90100

91101
% Option13: MAP probability curve
92-
out.option13 = listImagesStatSTEM('MAP probability','showMAPprob','outputMAP.MAPprob');
102+
out.option13 = listImagesStatSTEM('MAP probability','showMAPprob','outputMAP.MAPprob');

‎GUI/generalFunc/listStructNamesClass.m‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -23,4 +23,6 @@
2323
out(5,:) = {'atomcounting','atomCountStat'};
2424
out(6,:) = {'libcounting','atomCountLib'};
2525
out(7,:) = {'strainmapping','strainMapping'};
26-
out(8,:) = {'model3D','mod3D'};
26+
out(8,:) = {'model3D','mod3D'};
27+
out(9,:) = {'inputHMM','inputStatSTEM_HMM'};
28+
out(10,:) = {'outputHMM','outputStatSTEM_HMM'};

‎functions/@inputStatSTEM_HMM/analyseFactorialHMM.m‎

Lines changed: 113 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,113 @@
1+
function outputHMM = analyseFactorialHMM(inputHMM)
2+
3+
%%%% outputHMM = runEM(inputHMM);
4+
O = inputHMM.O;
5+
libraries = inputHMM.libraries;
6+
num_components = inputHMM.num_components;
7+
8+
% starting values
9+
ScalingHybrid = 1;
10+
Sigma = (max(O(:)) - min(O(:)))/2/(inputHMM.num_components);
11+
Init = ones(1,inputHMM.num_components)/(inputHMM.num_components);
12+
A = ones(inputHMM.num_components,inputHMM.num_components)/(inputHMM.num_components);
13+
14+
cyc = 1000;
15+
tol = 10^(-12);
16+
17+
H_viterbi = [];
18+
LL_viterbi = 0;
19+
coordinates_T = inputHMM.coordinates_T;
20+
outputHMM = outputStatSTEM_HMM(ScalingHybrid,Sigma,Init,A,H_viterbi,LL_viterbi,coordinates_T); %inputHMMect uit die klasse aanmaken eerst
21+
22+
loglikelihood = -Inf;
23+
24+
if isempty(inputHMM.GUI)
25+
outputHMM.GUI = inputHMM.GUI;
26+
end
27+
if ~isempty(inputHMM.GUI)
28+
inputHMM.waitbar.setValue(0)
29+
else
30+
h_bar = waitbar(0,'Fitting... ');
31+
end
32+
33+
34+
for cycle=1:cyc
35+
oldlikelihood = loglikelihood;
36+
37+
[gamma,xi,loglikelihood] = runEstep(O,libraries,num_components,ScalingHybrid,Sigma,Init,A);
38+
[A,Init,Sigma,ScalingHybrid] = runMstep(O,libraries,gamma,xi);
39+
40+
check = checkConvergence(oldlikelihood,loglikelihood,tol,cycle);
41+
if check == 1 || cycle == cyc
42+
outputHMM.maxLogLikelihood = loglikelihood;
43+
outputHMM.A = A;
44+
outputHMM.Sigma = Sigma;
45+
outputHMM.Init = Init;
46+
outputHMM.ScalingHybrid = ScalingHybrid;
47+
break;
48+
end
49+
50+
% Update waitbar
51+
if ~isempty(inputHMM.GUI)
52+
inputHMM.waitbar.setValue(cycle/cyc*100)
53+
else
54+
waitbar(cycle/cyc,h_bar,'Fitting ...')
55+
end
56+
57+
end
58+
% Update waitbar
59+
if ~isempty(inputHMM.GUI)
60+
inputHMM.waitbar.setValue(95)
61+
else
62+
delete(h_bar)
63+
end
64+
%%%% outputHMM = Viterbi(inputHMM,outputHMM);
65+
T = size(O,1);
66+
N = size(O,2);
67+
68+
delta = zeros(T,num_components,N);
69+
delta_all = zeros(T,num_components,num_components,N);
70+
arg = zeros(T,num_components,N);
71+
H = zeros(T,N);
72+
LL = zeros(1,N);
73+
74+
Gauss = calculateEmissionProbability(O,libraries,num_components,ScalingHybrid,Sigma);
75+
76+
% initialise best score and argument array
77+
delta(1,:,:) = Init.*reshape(Gauss(1,:,:),N,num_components)';
78+
arg(1,:,:) = 0;
79+
80+
% recursion
81+
for t=2:T
82+
for n = 1:N
83+
for i_g = 1:(num_components)
84+
delta_all(t,:,i_g,n) = delta(t-1,:,n).*A(:,i_g)';
85+
[delta(t,i_g,n),arg(t,i_g,n)] = max(reshape(delta_all(t,:,i_g,n),num_components,1));
86+
delta(t,i_g,n) = delta(t,i_g,n)*Gauss(t,n,i_g);
87+
end
88+
delta(t,:,n) = delta(t,:,n)/sum(delta(t,:,n)); %normalisatie owv computationele redenen!
89+
end
90+
end
91+
92+
% termination step at time T
93+
for n = 1:N
94+
[LL(n), H(T,n)] = max(delta(T,:,n)); %nog steeds LL door normalisatie owv computationele redenen?
95+
end
96+
LL = prod(LL);
97+
98+
% path backtracking
99+
for t = T-1:-1:1
100+
for n = 1:N
101+
H(t,n) = arg(t+1,H(t+1,n),n);
102+
end
103+
end
104+
105+
if libraries(1) == 0
106+
outputHMM.H_viterbi = H - 1;
107+
else
108+
outputHMM.H_viterbi = H;
109+
end
110+
111+
outputHMM.LL_viterbi = LL;
112+
113+
end

‎functions/@inputStatSTEM_HMM/inputStatSTEM_HMM.m‎

Lines changed: 93 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,93 @@
1+
classdef inputStatSTEM_HMM < StatSTEMfile
2+
% inputHMM - class to hold input parameters for fitting a hidden Markov
3+
% Model to a time series of ADF STEM images
4+
%
5+
6+
%--------------------------------------------------------------------------
7+
% This file is part of StatSTEM
8+
%
9+
% Copyright: 2019, EMAT, University of Antwerp
10+
% Author: A. De wael
11+
% License: Open Source under GPLv3
12+
% Contact: sandra.vanaert@uantwerpen.be
13+
%--------------------------------------------------------------------------
14+
15+
properties % input values
16+
libraries = []; % library of simulated cross-sections
17+
O = []; % observed sequence
18+
G = 0; % maximum number of atoms in a column
19+
num_components = 0;
20+
end
21+
22+
properties % additional input values for displaying results
23+
obs_T = []; % observations: all ADF STEM images
24+
model_T = []; % fitted models for all ADF STEM images
25+
% coordinates & dx is already in the StatSTEMfile class
26+
coordinates_T = []; % coordinates of all columns throughout all images
27+
end
28+
29+
methods
30+
outputHMM = analyseFactorialHMM(inputHMM);
31+
end
32+
33+
methods
34+
function obj = inputStatSTEM_HMM(O,libraries,G,obs_T,model_T,coordinates_T,dx)
35+
obj.O = O;
36+
obj.libraries = libraries;
37+
obj.G = G;
38+
if libraries(1) == 0
39+
obj.num_components = G + 1;
40+
else
41+
obj.num_components = G;
42+
end
43+
obj.obs_T = obs_T;
44+
obj.model_T = model_T;
45+
obj.coordinates_T = coordinates_T;
46+
obj.dx = dx;
47+
end
48+
49+
function obj = set.O(obj,val)
50+
% Check if input is not empty
51+
if isa(val,'double') && (isempty(val) || ( size(val,1)>=0 && size(val,2)>1 ) )
52+
obj.O = val;
53+
else
54+
error('StatSTEMfile: Observed sequence matrix O not properly defined')
55+
end
56+
end
57+
58+
function obj = set.libraries(obj,val)
59+
% Check if input is not empty
60+
if isa(val,'double') && (isempty(val) || ( size(val,1)>=0 && size(val,2)>=0 ) )
61+
if size(val,1) == 1 && size(val,2) > 1
62+
obj.libraries = val';
63+
elseif size(val,2) == 1 && size(val,1) > 1
64+
obj.libraries = val;
65+
end
66+
else
67+
error('StatSTEMfile: libraries not properly defined')
68+
end
69+
end
70+
71+
function obj = set.G(obj,val)
72+
% Check if input is not empty
73+
if isa(val,'double') && (isempty(val) || ( size(val,1)>=0 && size(val,2)==1 ) )
74+
obj.G = val;
75+
else
76+
error('StatSTEMfile: G not properly defined')
77+
end
78+
end
79+
80+
% function obj = set.obs_T(obj,val)
81+
% % Check if input is not empty
82+
% if isa(val,'double') && (isempty(val) || ( size(val,1)>=0 && size(val,2) >=0 && size(val,3) >=0 ) )
83+
% obj.obs_T = val;
84+
% else
85+
% error('StatSTEMfile: G not properly defined')
86+
% end
87+
% end
88+
89+
90+
end
91+
92+
93+
end

‎functions/@inputStatSTEM_HMM/showModel_T.m‎

Lines changed: 60 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,60 @@
1+
function showModel_T(obj)
2+
% showModel_T - Show the fitted model in the StatSTEM interface
3+
%
4+
% syntax: showModel_T(obj,t)
5+
% obj - inputStatSTEM object
6+
% t - frame number
7+
%
8+
9+
%--------------------------------------------------------------------------
10+
% This file is part of StatSTEM
11+
%
12+
% Copyright: 2018, EMAT, University of Antwerp
13+
% Author: K.H.W. van den Bos
14+
% License: Open Source under GPLv3
15+
% Contact: sandra.vanaert@uantwerpen.be
16+
%--------------------------------------------------------------------------
17+
18+
if isempty(obj.model_T)
19+
return
20+
end
21+
22+
T = size(obj.model_T,3);
23+
x_axis = (1:size(obj.model_T,2))*obj.dx;
24+
y_axis = (1:size(obj.model_T,1))*obj.dx;
25+
26+
tab = obj.ax.Parent.Parent.Parent;
27+
usr = get(tab,'Userdata');
28+
pan = obj.ax.Parent;
29+
imT = imagesc(x_axis,y_axis,obj.model_T(:,:,1));
30+
axis equal;
31+
axis off;
32+
colormap gray
33+
axis([x_axis(1),x_axis(end),y_axis(1),y_axis(end)]) % Set limits
34+
caxis([min(min(obj.obs_T(:))),max(max(obj.obs_T(:)))])
35+
36+
htext = 0.025;
37+
wtext = 0.05;
38+
hslider = 0.025;
39+
wslider = 0.2;
40+
posText = [0.5-wtext/2 0.1 wtext htext];
41+
posSlider = [posText(1)+wtext/2-wslider/2 posText(2)-htext*2 wslider hslider];
42+
TextT = uicontrol('Parent',pan,'style','text',...
43+
'units','normalized','position',posText);
44+
SliderT = uicontrol('Parent',pan,'style','slider','units','normalized','position',posSlider,...
45+
'min', 1, 'max', T,'Value',1);
46+
addlistener(SliderT, 'Value', 'PostSet', @callbackfn);
47+
function callbackfn(source,eventdata)
48+
t = round(get(eventdata.AffectedObject, 'Value'));
49+
50+
% Show observation
51+
imT.CData = obj.model_T(:,:,t);
52+
axis equal;
53+
axis off;
54+
colormap gray
55+
axis([x_axis(1),x_axis(end),y_axis(1),y_axis(end)]) % Set limits
56+
caxis([min(min(obj.model_T(:,:,t))),max(max(obj.model_T(:,:,t)))])
57+
58+
TextT.String = num2str(t);
59+
end
60+
end

0 commit comments

Comments
 (0)