improved variable naming | improved masking

Author: MSenden

code/matlab/IRM.m

@@ -25,8 +25,8 @@
     %
     % Input-referred model (IRM) mapping tool.
     %
-    % irm = IRM(params) creates an instance of the IRM class.
-    % params is a structure with 5 required fields
+    % irm = IRM(parameters) creates an instance of the IRM class.
+    % parameters is a structure with 5 required fields
     %   - f_sampling: sampling frequency (1/TR)
     %   - n_samples : number of samples (volumes)
     %   - n_rows    : number of rows (in-plane resolution)
@@ -55,7 +55,7 @@
     % (e.g. help IRM.mapping)
     %
     % typical workflow:
-    % 1. irm = IRM(params);
+    % 1. irm = IRM(parameters);
     % 2. irm.set_stimulus();
     % 3. irm.create_timecourse(FUN,xdata);
     % 4. results = irm.mapping(data);
@@ -87,24 +87,24 @@
 
     methods (Access = public)
 
-        function self = IRM(params,varargin)
+        function self = IRM(parameters,varargin)
             % constructor
             p = inputParser;
-            addRequired(p,'params',@isstruct);
+            addRequired(p,'parameters',@isstruct);
             addOptional(p,'hrf',[]);
-            p.parse(params,varargin{:});
+            p.parse(parameters,varargin{:});
 
             self.is = 'input-referred modeling tool';
 
             self.two_gamma = @(t) (6*t.^5.*exp(-t))./gamma(6)...
                 -1/6*(16*t.^15.*exp(-t))/gamma(16);
 
-            self.f_sampling = p.Results.params.f_sampling;
+            self.f_sampling = p.Results.parameters.f_sampling;
             self.p_sampling = 1/self.f_sampling;
-            self.n_samples = p.Results.params.n_samples;
-            self.n_rows = p.Results.params.n_rows;
-            self.n_cols = p.Results.params.n_cols;
-            self.n_slices = p.Results.params.n_slices;
+            self.n_samples = p.Results.parameters.n_samples;
+            self.n_rows = p.Results.parameters.n_rows;
+            self.n_cols = p.Results.parameters.n_cols;
+            self.n_slices = p.Results.parameters.n_slices;
             self.n_total = self.n_rows*self.n_cols*self.n_slices;
 
             if ~isempty(p.Results.hrf)
@@ -228,7 +228,7 @@ function create_timecourse(self,FUN,xdata)
             %
             % optional inputs are
             %  - threshold: minimum voxel intensity (default = 100.0)
-            % - mask      : binary mask for selecting voxels
+            %  - mask     : binary mask for selecting voxels
 
             text = 'mapping input-referred model...';
             fprintf('%s\n',text)
@@ -247,7 +247,7 @@ function create_timecourse(self,FUN,xdata)
             data = reshape(data(1:self.n_samples,:,:,:),...
                 self.n_samples,self.n_total);
             mean_signal = mean(data);
-            data = zscore(data);
+            
 
             if isempty(mask)
                 mask = mean_signal>=threshold;
@@ -256,53 +256,51 @@ function create_timecourse(self,FUN,xdata)
             voxel_index = find(mask);
             n_voxels = numel(voxel_index);
 
-            mag_d = sqrt(sum(data(:,mask).^2));
+            data = zscore(data(:,mask));
+            mag_d = sqrt(sum(data.^2));
 
             results.R = zeros(self.n_total,1);
             results.P = zeros(self.n_total,self.n_predictors);
 
             if size(self.hrf,2)==1
-                hrf_fft = fft(repmat([self.hrf;...
-                    zeros(self.n_samples-self.l_hrf,1)],...
-                    [1,self.n_points]));
+                hrf_fft = fft([self.hrf;...
+                    zeros(self.n_samples-self.l_hrf,1)]);
                 tc = zscore(ifft(self.tc_fft.*hrf_fft))';
 
                 mag_tc = sqrt(sum(tc.^2,2));
                 for m=1:n_voxels
                     v = voxel_index(m);
 
-                    CS = (tc*data(:,v))./...
-                        (mag_tc*mag_d(v));
+                    CS = (tc*data(:,m))./...
+                        (mag_tc*mag_d(m));
                     id = isinf(CS) | isnan(CS);
                     CS(id) = 0;
                     [results.R(v),j] = max(CS);
                     for p=1:self.n_predictors
                         results.P(v,p) = self.xdata{p}(self.idx(j,p));
                     end
 
-                    waitbar(v/self.n_total,wb)
+                    waitbar(m/n_voxels,wb)
                 end
             else
-                hrf_fft_all = fft([self.hrf;...
-                    zeros(self.n_samples-self.l_hrf,self.n_total)]);
+                hrf_fft_all = fft([self.hrf(:,mask);...
+                    zeros(self.n_samples-self.l_hrf,n_voxels)]);
                 for m=1:n_voxels
                     v = voxel_index(m);
 
-                    hrf_fft = repmat(hrf_fft_all(:,v),...
-                        [1,self.n_points]);
-                    tc = zscore(ifft(self.tc_fft.*hrf_fft))';
+                    tc = zscore(ifft(self.tc_fft.*hrf_fft_all(:,m)))';
                     mag_tc = sqrt(sum(tc.^2,2));
 
-                    CS = (tc*data(:,v))./...
-                        (mag_tc*mag_d(v));
+                    CS = (tc*data(:,m))./...
+                        (mag_tc*mag_d(m));
                     id = isinf(CS) | isnan(CS);
                     CS(id) = 0;
                     [results.R(v),j] = max(CS);
                     for p=1:self.n_predictors
                         results.P(v,p) = self.xdata(self.idx(j,p),p);
                     end
 
-                    waitbar(v/self.n_total,wb)
+                    waitbar(m/n_voxels,wb)
                 end
             end
             results.R = reshape(results.R,self.n_rows,self.n_cols,self.n_slices);