Newer
Older
classdef Interface < handle
properties
T %temps (ou numéro de frame)
Y %observation : volume du ventricule gauche. Une observation (frame) par ligne de Y
N
meilleure_sequence
meilleure_sequence_models
loglikelihood
%% morceaux courants (cad non finis) A, B1, B2 et C
dynamics_models
first_frame_of_models
sigma_estimation
sigma_test
proba_morceaux % probabilité courante pour chaque morceau, avec son meilleur path Viterbi
morceaux_finis % pour chaque morceau courant : liste des morceaux finis du path. Cette liste contient une copie du morceau dans l'état où il était quand il a été terminé
proba_morceaux_finis % pour chaque morceau courant : liste des probabilités associées aux morceaux finis. Elle combine la proba de chaque observation du morceau et la proba de transition du morceau au suivant
model_indexes % pour chaque morceau courant : liste des indices des morceaux du path Viterbi pour chaque frame
scores_finis % pour chaque morceau courant : liste des proba d'observation pour chaque frame (uniquement les frames des morceaux finis, car le morceau courant n'est pas encore stable et les x peuvent encore changer)
end
methods
function obj = Interface(interpolant_filenames, display_screen)
obj.interpolant_filenames = interpolant_filenames;
obj.transition = [ 0.5 0.5 0 0; % départ de A
0 1/3 1/3 1/3; % départ de B1
0 0 0.5 0.5; % départ de B2
0 0 0 1 % départ de C
];
obj.display_screen = display_screen;
end
function process_new_sequence(obj, filename, initial_alphas, sigma_estimation, sigma_test, threshold_dynamics)
obj.threshold_dynamics = threshold_dynamics;
obj.dynamics_models = {NaN, NaN, NaN, NaN}; % modèles courants A B1 B2 et C
obj.first_frame_of_models = [1 NaN NaN NaN]; % frame où chaque modèle commence
if obj.display_screen == -1
obj.dynamics_models{1} = Dynamics_Model(obj.interpolant_filenames{1}, obj.display_screen);
else
obj.dynamics_models{1} = Dynamics_Model(obj.interpolant_filenames{1}, obj.display_screen+1);
end
obj.proba_morceaux = [1 0 0 0]; % proba des modèles A, B1, B2, et C
obj.T = [];
obj.Y = [];
obj.N = [];
obj.scores_finis = {[], [], [], []};
obj.model_indexes = {[], [], [], []};
obj.initial_alphas = initial_alphas;
obj.sigma_estimation = sigma_estimation;
obj.sigma_test = sigma_test;
obj.morceaux_finis = {[], [], [], []};
obj.proba_morceaux_finis = {[], [], [], []};
obj.read_data(filename);
obj.N = size(obj.T,1);
obj.dynamics_models{1}.initialise(0, initial_alphas(1), sigma_estimation, sigma_test);
%ligne 14
for frame=1:obj.N
obj.process_frame(frame);
if obj.display_screen ~= -1
for model=1:4
if obj.proba_morceaux(model) == 0 || isnan(obj.proba_morceaux(model))
continue
end
obj.dynamics_models{model}.plot_estimated_X();
end
end
end
%% choix de la meilleure sequence (peut être pour plus tard: toujours choisir la séquence qui se termine sur le modèle C ?)
[max_proba, indice] = max(obj.proba_morceaux);
obj.meilleure_sequence = obj.model_indexes{indice}
obj.meilleure_sequence_models = [obj.morceaux_finis{indice} obj.dynamics_models{indice}];
%% calcul les scores pour les frames du dernier morceau
new_scores = obj.meilleure_sequence_models(end).compute_loglikelihood_dynamics([1:obj.N-obj.first_frame_of_models(indice)+1]);
obj.loglikelihood = [obj.scores_finis{indice} new_scores];
%% à faire : afficher la séquence des modèles et la séquence de scores comme dans plot
close all
for i=1:length(obj.meilleure_sequence_models)
obj.meilleure_sequence_models(i).plot_estimated_X();
end
%%% save plot
end
function read_data(obj, filename)
obj.Y =load(filename);
%obj.T = [1:75]';
%obj.Y = obj.Y / max(obj.Y);
%obj.Y = repmat(obj.Y, 3, 1);
obj.T = [1:25]';
obj.Y = obj.Y / max(obj.Y);
end
function process_frame(obj, frame_number)
frame_number
Yt = obj.Y(frame_number,:);
Tt = obj.T(frame_number);
%% tableaux temporaires dans lesquels on prépare la prochaine version des tableaux de la classe sans les modifier avant d'avoir terminé de traiter la frame pour tous les morceaux
new_proba_morceaux = [NaN NaN NaN NaN];
new_morceaux_finis = {[], [], [], []};
new_proba_morceaux_finis = {[], [], [], []};
new_models = {NaN, NaN, NaN, NaN};
new_first_frame_of_models = [NaN NaN NaN NaN];
new_scores_finis = {[], [], [], []};
new_model_indexes = {[], [], [], []};
%% boucle sur les 4 modèles possibles
for morceau = 1:4
% on va regarder toutes les branches qui arrivent sur le modèle
%% tableaux temporaires pour le morceau courant et toutes les branches possibles
modele_branche = {NaN, NaN, NaN, NaN}; % copie des modèles pour leur montrer la nouvelle frame, sans modifier ceux de obj.dynamics_models
first_frame_of_model_branche = [NaN NaN NaN NaN]; % pour les copies des modèles
proba_obs_branche = [NaN NaN NaN NaN]; % proba d'observation donnée par les copies des modèles ayant vu la frame courante
nouveau_morceau_fini_branche = {NaN, NaN, NaN, NaN}; % si une branche implique un nouveau morceau fini, on le sauvegarde ici en attendant Viterbi et le choix de la branche (cad si on le garde ou pas)
proba_branche = [NaN NaN NaN NaN]; % proba totale de chaque branche
proba_morceau_prec = [NaN NaN NaN NaN]; % proba totale des frames vues par le morceau qui était courant à la frame précédente (donc d'indice=branche). Ces proba sont combinées avec les probas de transition au sein du morceau
new_scores_branche = {[] [] [] []}; % score des frames vues par le morceau qui était courant à la frame précédente
new_models{morceau} = copy(obj.dynamics_models{morceau});
new_first_frame_of_models(morceau) = obj.first_frame_of_models(morceau);
new_models{morceau}.process_new_frame(Yt, Tt);
proba_obs = new_models{morceau}.compute_loglikelihood_dynamics(1);
new_proba_morceaux(morceau) = log(obj.proba_morceaux(morceau)) + proba_obs;
new_morceaux_finis{morceau} = [];
new_model_indexes{morceau} = [morceau];
else
%% on calcule les probas de toutes les branches qui arrivent sur le morceau courant
for branche = 1:4
trans = obj.transition(branche, morceau);
% on commence par ignorer les branches impossibles
if trans == 0
continue
end
if obj.proba_morceaux(branche) == 0 || isnan(obj.proba_morceaux(branche))
continue
end
if branche == morceau % cas sans changement de modèle, donc le modèle est déjà initialisé
modele_branche{branche} = copy(obj.dynamics_models{branche});
first_frame_of_model_branche(branche) = obj.first_frame_of_models(morceau);
else % cas avec changement de modèle, donc il faut créer et initialiser le modèle
%% on crée le modèle et on l'initialise
if obj.display_screen == -1
modele_branche{branche} = Dynamics_Model(obj.interpolant_filenames{morceau}, obj.display_screen);
else
modele_branche{branche} = Dynamics_Model(obj.interpolant_filenames{morceau}, obj.display_screen+morceau);
end
modele_branche{branche}.initialise(0, obj.initial_alphas(morceau), obj.sigma_estimation, obj.sigma_test);
%modele_branche{branche}.initialise(0, obj.dynamics_models{branche}.current_alpha, obj.sigma_estimation, obj.sigma_test);
first_frame_of_model_branche(branche) = frame_number;
end
%% on donne la frame au modèle, et on récupère la proba d'observation (dans la fenêtre)
modele_branche{branche}.process_new_frame(Yt, Tt);
%print branche et morceau
BRANCHE=branche
MORCEAU=morceau
proba_obs_branche(branche) = modele_branche{branche}.compute_loglikelihood_dynamics(frame_number-first_frame_of_model_branche(branche)+1);
%% on garde en mémoire les morceaux qui sont finis, pour l'affichage final des X
%%%% %%%%%%%%%%%%%%%% a faire : mettre à jour les x pour bien les étaler
if branche == morceau
nouveau_morceau_fini_branche{branche} = 0;
nouveau_morceau_fini_branche{branche} = copy(obj.dynamics_models{branche}); % ce modèle n'a pas encore vu la frame courante, on peut donc le sauvegader dans son état à la frame précédente
%%%% %%%%%% mettre à jour la copie ici
%% on récupère les scores des frames vues par le modèle
if branche == morceau
%% on récupère les scores pour les frames précédentes, qui ont pu bouger avec le argmax
proba_morceau_prec(branche) = sum(modele_branche{branche}.compute_loglikelihood_dynamics([1:frame_number-first_frame_of_model_branche(branche)]));
proba_morceau_prec(branche) = proba_morceau_prec(branche) + log(trans) * (frame_number-first_frame_of_model_branche(branche)-1);
else
%% on récupère les scores pour les frames précédentes, qui ont pu bouger avec le argmax
new_scores_branche{branche} = nouveau_morceau_fini_branche{branche}.compute_loglikelihood_dynamics([1:frame_number-obj.first_frame_of_models(branche)]);
proba_morceau_prec(branche) = sum(new_scores_branche{branche});
trans_morceau_prec = obj.transition(branche, branche);
proba_morceau_prec(branche) = proba_morceau_prec(branche) + log(trans_morceau_prec) * (frame_number-obj.first_frame_of_models(branche)-1);
proba_branche(branche) = sum(obj.proba_morceaux_finis{branche}) + proba_morceau_prec(branche) + log(trans) + proba_obs_branche(branche);
end
%% on choisi la meilleure branche
[max_val, indice] = max(proba_branche);
if isnan(max_val)
continue
end
new_proba_morceaux(morceau) = max_val;
new_models{morceau} = modele_branche{indice};
new_first_frame_of_models(morceau) = first_frame_of_model_branche(indice);
if nouveau_morceau_fini_branche{indice} ~= 0
new_morceaux_finis{morceau} = [obj.morceaux_finis{indice} nouveau_morceau_fini_branche{indice}];
trans = obj.transition(indice, morceau);
new_proba_morceaux_finis{morceau} = [obj.proba_morceaux_finis{indice} proba_morceau_prec(indice) + log(trans)];
new_scores_finis{morceau} = [obj.scores_finis{indice} new_scores_branche{indice}];
else
new_morceaux_finis{morceau} = obj.morceaux_finis{indice};
new_proba_morceaux_finis{morceau} = obj.proba_morceaux_finis{indice};
new_scores_finis{morceau} = obj.scores_finis{indice};
new_model_indexes{morceau} = [obj.model_indexes{indice} morceau];
obj.proba_morceaux = new_proba_morceaux;
obj.dynamics_models = new_models;
obj.first_frame_of_models = new_first_frame_of_models;
obj.morceaux_finis = new_morceaux_finis;
obj.proba_morceaux_finis = new_proba_morceaux_finis;
obj.model_indexes = new_model_indexes;
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
%% réflexions pour plus tard
%% gérer : alpha différents pour les différents modèles (alpha_ini + marges)
%% fin des modèles : duppliquer les dernières colonnes pour que x ait le temps de se fixer ? (avec pénalité lors de la comparaison des modèles)
%% ou alors : on n'attend pas que les x soient stables
%% initialiser le modèle concurrent à chaque fois (tant qu'on n'a pas basculé dessus)
end
function plot(obj) %% cette fonction est à modifier car maintenant nous avons plusieurs modèles / morceaux
obj.dynamics_model.plot_estimated_X();
figure(obj.display_screen);
n = length(obj.loglikelihood_dynamics_window);
subplot(4,1,1)
hold off
plot(obj.T(1:n),obj.Y(1:n,1))
title('Pose representation')
xlabel('frame number')
ylabel('1st dimension')
XL = xlim;
subplot(4,1,2)
hold off
plot(obj.T(1:n),obj.loglikelihood_pose)
hold on
indexes_p = (obj.loglikelihood_pose < obj.threshold_pose);
plot(obj.T(indexes_p),obj.loglikelihood_pose(indexes_p),'.r')
title('Likelihood of the pose to be normal')
xlabel('Frame number')
ylabel('Log likelihood')
xlim(XL)
YL = ylim;
YL1 = max(YL(1),-5);
ylim([YL1 YL(2)])
hold on
subplot(4,1,3)
hold off
plot(obj.T(1:n),obj.loglikelihood_dynamics_window)
hold on
indexes_w = (obj.loglikelihood_dynamics_window < obj.threshold_dynamics);
plot(obj.T(indexes_w),obj.loglikelihood_dynamics_window(indexes_w),'.r')
title('Likelihood of the dynamics to be normal, computed in the dynamic window')
xlabel('Frame number')
ylabel('Log likelihood')
xlim(XL)
YL = ylim;
YL1 = max(YL(1),-20);
ylim([YL1 YL(2)])
hold on
subplot(4,1,4)
hold off
n2 = length(obj.loglikelihood_dynamics);
plot(obj.T(1:n2),obj.loglikelihood_dynamics)
hold on
indexes_d = (obj.loglikelihood_dynamics < obj.threshold_dynamics);
plot(obj.T(indexes_d),obj.loglikelihood_dynamics(indexes_d),'.r')
title('Likelihood of the dynamics to be normal, computed for each converged frame')
xlabel('Frame number')
ylabel('Log likelihood')
xlim(XL)
YL = ylim;
YL1 = max(YL(1),-20);
YL2 = min(YL(2),10);
ylim([YL1 YL2])
hold on
subplot(4,1,1)
hold on
plot(obj.T(indexes_p),obj.Y(indexes_p,1),'.r')
%plot(obj.T(indexes_w),obj.Y(indexes_w,1),'.r')
plot(obj.T(indexes_d),obj.Y(indexes_d,1),'.r')
%indexes = indexes_p == 0 & indexes_w == 0 & [indexes_d ones(1,n-n2)] == 0;
indexes = indexes_p == 0 & [indexes_d ones(1,n-n2)] == 0;
plot(obj.T(indexes),obj.Y(indexes,1),'.g')
end
end
end