Merge pull request #142 from sampsapursiainen/main_development_branch

Zeffiro developer: A regular push adding the changes made in the curr…

+zeffiro/+plugins/+Kalman/+m/RTS_smoother.m

@@ -1,34 +1,3 @@
 function [P_s_store, m_s_store, G_store] = RTS_smoother(P_store, z_inverse, A, Q, number_of_frames)
-P_s_store = cell(0);
-m_s_store = cell(0);
-G_store = cell(0);
-h = zef_waitbar(0,'Smoothing');
-for f_ind = number_of_frames:-1:1
-    zef_waitbar(1 - f_ind/number_of_frames,h, ['Smoothing ' int2str(number_of_frames -f_ind) ' of ' int2str(number_of_frames) '.']);
-
-    P = P_store{f_ind};
-    m = z_inverse{f_ind};
-    % if A is Identity
-    if (isdiag(A) && all(diag(A) - 1) < eps)
-        P_ = P + Q;
-        m_ = m;
-        G =  P / P_;
-    else
-        P_ = A * P * A' + Q;
-        m_ = A * m;
-        G =  (P * A') / P_;
-    end
-    if f_ind == number_of_frames
-        m_s = m;
-        P_s = P;
-    else
-        m_s = m + G * (m_s - m_);
-        P_s = P + G * (P_s - P_) * G';
-    end
-    P_s_store{f_ind} = P_s;
-    G_store{f_ind} = G;
-    m_s_store{f_ind} = m_s;
-end
-
-close(h);
+[P_s_store, m_s_store, G_store] = RTS_smoother(P_store, z_inverse, A, Q, number_of_frames);
 end

+zeffiro/+plugins/+Kalman/+m/zef_KF.m

@@ -1,94 +1,3 @@
-function [zef] = zef_KF(zef, q_value)
-% Optimal q_value as parameter
-%% Initial parameters
-snr_val = zef.inv_snr;
-pm_val = zef.inv_prior_over_measurement_db;
-amplitude_db = zef.inv_amplitude_db;
-pm_val = pm_val - amplitude_db;
-std_lhood = 10^(-snr_val/20);
-sampling_freq = zef.inv_sampling_frequency;
-high_pass = zef.inv_low_cut_frequency;
-low_pass = zef.inv_high_cut_frequency;
-number_of_frames = zef.number_of_frames;
-source_direction_mode = zef.source_direction_mode; 
-source_directions = zef.source_directions;
-source_positions = zef.source_positions;
-time_step = zef.inv_time_3;
-
-%% Reconstruction identifiers
-reconstruction_information.tag = 'Kalman';
-reconstruction_information.inv_time_1 = zef.inv_time_1;
-reconstruction_information.inv_time_2 = zef.inv_time_2;
-reconstruction_information.inv_time_3 = zef.inv_time_3;
-reconstruction_information.sampling_freq = zef.inv_sampling_frequency;
-reconstruction_information.low_pass = zef.inv_high_cut_frequency;
-reconstruction_information.high_pass = zef.inv_low_cut_frequency;
-reconstruction_information.number_of_frames = zef.number_of_frames;
-reconstruction_information.source_direction_mode = zef.source_direction_mode;
-reconstruction_information.source_directions = zef.source_directions;
-reconstruction_information.snr_val = zef.inv_snr;
-reconstruction_information.pm_val = zef.inv_prior_over_measurement_db;
-
-%%
-[L,n_interp, procFile] = zef_processLeadfields(zef);
-
-%get ellipse filteres full measurement data. f_data: "sensors" x "time points"
-[f_data] = zef_getFilteredData; 
-timeSteps = arrayfun(@(x) zef_getTimeStep(f_data, x, zef), 1:number_of_frames, 'UniformOutput', false);
-
-z_inverse_results = cell(0);
-%% CALCULATION STARTS HERE
-% m_0 = prior mean
-m = zeros(size(L,2), 1);
-
-[theta0] = zef_find_gaussian_prior(snr_val-pm_val,L,size(L,2),zef.normalize_data,0);
-
-% Transition matrix is Identity matrix
-P = eye(size(L,2)) * theta0;
-
-A = eye(size(L,2));
-
-% If q_value given in the function call
-if nargin > 1
-    Q = q_value*eye(size(L,2));
-else
-    zef_init_gaussian_prior_options;
-    evolution_prior_db = zef.inv_evolution_prior;
-    q_value = find_evolution_prior(L, theta0, time_step, evolution_prior_db);
-    Q = q_value*eye(size(L,2));
-end
-reconstruction_information.Q = q_value;
-
-% std_lhood
-R = std_lhood^2 * eye(size(L,1));
-
-Q_Store = cell(0);
-
-%% KALMAN FILTER
-filter_type = zef.KF.filter_type.Value;
-smoothing = zef.kf_smoothing;
-if filter_type == '1'
-    [P_store, z_inverse] = zeffiro.plugins.Kalman.m.kalman_filter(m,P,A,Q,L,R,timeSteps, number_of_frames, smoothing);
-elseif filter_type == '2'
-    n_ensembles = str2double(zef.KF.number_of_ensembles.Value);
-    z_inverse = zeffiro.plugins.Kalman.m.EnKF(m,A,P,Q,L,R,timeSteps,number_of_frames, n_ensembles);
-elseif filter_type == '3'
-    [P_store, z_inverse] = zeffiro.plugins.Kalman.m.kalman_filter_sLORETA(m,P,A,Q,L,R,timeSteps, number_of_frames, smoothing);
-end
-
-%% RTS SMOOTHING
-
-if (smoothing == 2)
-[~, m_s_store, ~] = zeffiro.plugins.Kalman.m.RTS_smoother(P_store, z_inverse, A, Q, number_of_frames);
-z_inverse = m_s_store;
-end
-
-%% POSTPROCESSING
-[z] = zef_postProcessInverse(z_inverse, procFile);
-%normalize the reconstruction so that the highest value is equal to 1
-[z] = zef_normalizeInverseReconstruction(z);
-%% CALCULATION ENDS HERE
-zef.reconstruction_information = reconstruction_information;
-zef.reconstruction = z;
-
+function [z, reconstruction_information] = zef_KF(q_value)
+[z, reconstruction_information] = zef_KF(q_value);
 end

m/lead_field/zef_eeg_lead_field.m

@@ -1,5 +1,9 @@
 function zef = zef_eeg_lead_field(zef)
 
+if nargin == 0
+zef = evalin('base','zef');
+end
+
 warning('off');
 zef.lead_field_type = 1;
 zef.imaging_method = 1;
@@ -13,4 +17,8 @@
 end
 warning('on');
 
+if nargout == 0
+assignin('base','zef',zef);
+end
+
 end

m/lead_field/zef_source_interpolation.m

@@ -1,7 +1,11 @@
 
 %Copyright © 2018- Sampsa Pursiainen & ZI Development Team
 %See: https://github.com/sampsapursiainen/zeffiro_interface
-function [source_interpolation_ind] = zef_source_interpolation(zef)
+function zef = zef_source_interpolation(zef)
+
+if nargin == 0
+zef = evalin('base','zef');
+end
 
     if evalin('base','isequal(size(zef.L,2),size(zef.source_directions,1))')
     evalin('base','zef.source_directions=zef.source_directions(find(not(isnan(sum(abs(zef.L),1)))),:);');
@@ -121,6 +125,8 @@
 MdlKDT = KDTreeSearcher(center_points);
 source_interpolation_ind{3} = knnsearch(MdlKDT,source_positions);
 
+zef.source_interpolation_ind = source_interpolation_ind;
+
 if nargout == 0
 assignin('base','zef', zef);
 end

m/zef_close_all.m

@@ -1,3 +1,11 @@
+function  zef_close_all(zef)
+
+if nargin == 0
+    if evalin('base','exist(''zef'',''var'')')
+    zef = evalin('base','zef');
+    end
+end
+
 zef.h_window_aux = findall(groot,'-regexp','Name','ZEFFIRO Interface*');
 %zef.h_window_aux = findall(groot,'-property','ZefTool','-or','-property','ZefFig');
 set(zef.h_window_aux,'DeleteFcn','');
@@ -21,6 +29,7 @@
 warning('on','MATLAB:rmpath:DirNotFound');    
 end
 
-clear zef zef_data zef_i zef_j zef_k
+evalin('base','clear zef zef_data zef_i zef_j zef_k;');
 
+end
 

scripts/KalmanDemo.m

@@ -24,4 +24,6 @@ function KalmanDemo()
 % Save project
 project_struct = zef_save(project_struct,'example_project.mat','data/');
 
+zef_close_all;
+
 end

scripts/zef_rec_max_diff.m

@@ -1,4 +1,4 @@
-function [dist_vec,angle_vec] = zef_rec_max_diff(zef,method_name)
+function [dist_vec,angle_vec,mag_vec] = zef_rec_max_diff(zef,method_name)
 
 dist_vec = zeros(3*size(zef.source_positions,1),1);
 angle_vec = zeros(3*size(zef.source_positions,1),1);
@@ -36,7 +36,7 @@
       for j = 1 : 3
 
    z_norm = sqrt(sum(reshape(z{3*(i-1) + j},3,size(zef.source_positions,1)).^2));
-   [~, I] = max(z_norm);
+   [mag_val, I] = max(z_norm);
 
    dir_vec_rec = z{3*(i-1) + j}(3*(I-1)+1:3*(I-1)+3);
    dir_vec_rec = dir_vec_rec/norm(dir_vec_rec,2);
@@ -45,6 +45,7 @@
 
    dist_vec(3*(i-1) + j) = (1/sqrt(3))*sqrt(sum((zef.source_positions(i,:) - zef.source_positions(I,:)).^2,2));
    angle_vec(3*(i-1) + j) = acosd(dot(dir_vec_rec,dir_vec_source));
+   mag_vec(3*(i-1) + j) = mag_val;
 
 
       end