add outlier functions to bypass matlab's isoutlier; adapt to r2017b

tICA/scripts/ComputeTICAFeatures.m

@@ -317,7 +317,7 @@
 %Find strongest individual subjects for components
 [~, I]=max(TCSVARS);
 for i=1:size(TCSVARS,2)
-    %if ~isfile([OutputFolder '/tICA' num2str(i,'%02.f') '_SS.sdseries.nii']) && ~isfile([OutputFolder '/All' num2str(i,'%02.f') '_SS.sdseries.nii'])
+    if ~isfile([OutputFolder '/tICA' num2str(i,'%02.f') '_SS.sdseries.nii']) && ~isfile([OutputFolder '/All' num2str(i,'%02.f') '_SS.sdseries.nii'])
         %Subjlist{I(i)}
         SubjFolderlist=[StudyFolder '/' Subjlist{I(i)}];
         tICAMaps_sub=ciftiopen([SubjFolderlist '/MNINonLinear/fsaverage_LR' LowResMesh 'k/' Subjlist{I(i)} '.' tICAFeaturesProcString '_SR' RegString '.' LowResMesh 'k_fs_LR.dscalar.nii'],'wb_command');      
@@ -355,14 +355,10 @@
         temp=tICAtcs_SS;
         temp.cdata=zeros(size(rtICASingleComponentSpaceXTime,1),size(tICAtcs_SS.cdata,2),'single');
         temp.cdata(:,1:size(rtICASingleComponentSpaceXTime,2))=squeeze(rtICASingleComponentSpaceXTime(:,:,i))./repmat(std(sum(rtICASingleComponentSpaceXTime,3),[],2),1,size(rtICASingleComponentSpaceXTime,2));
-        if ~isfile([OutputFolder '/tICA' num2str(i,'%02.f') '_SS.sdseries.nii'])
-            ciftisavereset(temp,[OutputFolder '/tICA' num2str(i,'%02.f') '_SS.sdseries.nii'],'wb_command');
-        end
+        ciftisavereset(temp,[OutputFolder '/tICA' num2str(i,'%02.f') '_SS.sdseries.nii'],'wb_command');
         temp.cdata(:,1:size(rtICASingleComponentSpaceXTime,2))=squeeze(sum(rtICASingleComponentSpaceXTime(:,:,:),3))./repmat(std(sum(rtICASingleComponentSpaceXTime,3),[],2),1,size(rtICASingleComponentSpaceXTime,2));
-        if ~isfile([OutputFolder '/All' num2str(i,'%02.f') '_SS.sdseries.nii'])
-            ciftisavereset(temp,[OutputFolder '/All' num2str(i,'%02.f') '_SS.sdseries.nii'],'wb_command');
-        end
-    %end
+        ciftisavereset(temp,[OutputFolder '/All' num2str(i,'%02.f') '_SS.sdseries.nii'],'wb_command');
+    end
 end
 if ~isfile([OutputFolder '/tICA_Maps_' num2str(tICAdim) '_' nonlinear '_SS.dscalar.nii'])
     ciftisave(tICAMaps_SS,[OutputFolder '/tICA_Maps_' num2str(tICAdim) '_tanhF_SS.dscalar.nii'],wbcommand);

tICA/scripts/feature_helpers/findFivePercentageOutliers.m

@@ -0,0 +1,12 @@
+function ret = findFivePercentageOutliers(flatdata)
+    %isoutlier(..., 'percentiles', [5 95]) written to treat all values at once
+    numelems = numel(flatdata);
+    ret = false(size(flatdata));
+    if numelems < 1
+        return
+    end
+    numexclude = floor((numelems + 9) / 20); %in R2021a, isoutlier excludes 0 for 10, 2 for 11...
+    [~, perm] = sort(flatdata(:));
+    ret([1:numexclude (end - numexclude + 1):end]) = true;
+    ret(perm) = ret; %un-sort the exclusion mask
+end

tICA/scripts/feature_helpers/findOutliers.m

@@ -0,0 +1,7 @@
+function ret = findOutliers(data)
+% must be 1d array
+% outlier detection: identify data that do not pass the MAD-based test, while ensuring that it never detects more than 5% of the highest or lowest values.
+madOutliers=findScaledMADOutliers(data);
+perOutliers=findFivePercentageOutliers(data);
+ret=and(madOutliers, perOutliers);
+end

tICA/scripts/feature_helpers/findScaledMADOutliers.m

@@ -0,0 +1,15 @@
+function ret = findScaledMADOutliers(data)
+% Determine if an array contains values that are outside of three scaled MAD
+    numElems = numel(data);
+    if numElems < 1
+        return
+    end
+    deviations = data - median(data);
+    absDeviations = abs(deviations);
+    mad = median(absDeviations);
+    % Scale the MAD by -1/(sqrt(2)*erfcinv(3/2))
+    % from https://www.mathworks.com/help/matlab/ref/isoutlier.html#bvolffm
+    scaledMAD = 1.4826 * mad;
+    threshold = 3 * scaledMAD;
+    ret = absDeviations > threshold;
+end

tICA/scripts/feature_helpers/std_outlier.m

@@ -7,7 +7,7 @@
     if any(isnan(data_reshape))
         error('found NaN in std_outlier');
     end
-    isoutlier_stat = myoutliers(data_reshape); %TODO: test me
+    isoutlier_stat = findOutliers(data_reshape); %TODO: test me
     if isempty(isoutlier_stat)
         %outlier_stat=1; %original
         outlier_stat=0;
@@ -18,17 +18,3 @@
         outlier_stat=1-std_exclude_outlier/std_all;
     end
 end
-
-%isoutlier(..., 'percentiles', [5 95]) written to treat all values at once
-function ret = myoutliers(flatdata)
-    numelems = numel(flatdata);
-    ret = false(size(flatdata));
-    if numelems < 1
-        return
-    end
-    numexclude = floor((numelems + 9) / 20); %in R2021a, isoutlier excludes 0 for 10, 2 for 11...
-    [~, perm] = sort(flatdata(:));
-    ret([1:numexclude (end - numexclude + 1):end]) = true;
-    ret(perm) = ret; %un-sort the exclusion mask
-end
-

tICA/scripts/feature_helpers/sum_outlier.m

@@ -1,8 +1,9 @@
 function output = sum_outlier(T)
 %Afill = filloutliers(T,'spline');
 %Afill = filloutliers(T,'next');
-Afill = filloutliers(T,'spline','percentiles', [5,95]);
+outliers=findOutliers(T);
 
+Afill = filloutliers(T,'spline','OutlierLocations', outliers);
 Afill(find(isnan(Afill)))=T(find(isnan(Afill)));
 % original
 output = sqrt(sum((T-Afill).^2));