Revert "Revert "Merge branch 'master' of github.com:canlab/CanlabCore""

This reverts commit a3f4b62.

CanlabCore/@atlas/downsample_parcellation.m

@@ -0,0 +1,134 @@
+function new_atlas_obj = downsample_parcellation(obj, varargin)
+% An atlas_obj has multiple label fields meant to label nested parcellations
+% at different levels of granularity. labels is the default which must
+% contain unique entries and correspond to the number of unique indices
+% in the (voxelwise) parcellation index map, but labels_2 ... labels_5
+% may contain non-unique entries corresponding to coarser parcellations.
+% This function remaps indices and probability maps to correspond to one
+% of these coarser parcellations. It is a lossy operations and all labels
+% of finer parcellation than the selected level will be discarded.
+%
+% Input ::
+%     obj - an atlas_object. Must have at least one of labels_2 ...
+%                 labels_5 defined
+%     varargin - either a vector of new labels to use, a string or leave empty. 
+%                  If you provide a vector, its length should equal number of
+%                  regions in original atlas, with redundant labels for regions
+%                  you want to combine in the downsampled atlas.
+%                  If you provide a string, it should be one of 
+%                  {'labels_2', 'labels_3', 'labels_4', 'labels_5'} which will
+%                  then be used as the downsampling vector.
+%                  If you don't provide anything the contents of labels_2 is used
+%                  by default.
+%
+% Output ::
+%     atlas_obj - a new atlas object, reindexed according to labelfield
+
+fprintf('Downsampling %s parcels\n', obj.atlas_name);
+
+if isempty(varargin)
+    labelfield = 'labels_2';
+elseif ischar(varargin{1})
+    labelfield=varargin{1};
+elseif isvector(varargin{1})
+    if length(varargin{1}) ~= num_regions(obj)
+        error('New labels has length %d which does not match input atlas parcel count (%d)',length(varargin{1}), num_regions(obj));
+    end
+    labelfield='labels';
+    obj.labels = varargin{1};
+else
+    error('Unrecognized input datatype');
+end
+
+% input check
+valid_lbls = {'labels', 'labels_2', 'labels_3', 'labels_4', 'labels_5'};
+if ~ismember(labelfield, valid_lbls)
+    error('labelfield must be one of {''labels_2'', ''labels_3'', ''labels_4'', ''labels_5''})');
+end
+
+% check that labels are ordered with finer labels first
+n_unique_lbls = zeros(1,length(valid_lbls));
+n_lbls = zeros(1,length(valid_lbls));
+for i = 1:length(valid_lbls)
+    lbl = valid_lbls{i};
+    n_unique_lbls(i) = length(unique(obj.(lbl)));
+    n_lbls(i) = length(obj.(lbl));
+end
+
+new_ind = find(strcmp(valid_lbls,labelfield));
+if n_unique_lbls(new_ind) > n_unique_lbls(1)
+    error('labelfield %s has %d labels, which is finer grained than labels which has %d unique labels',labelfield, n_unique_lbls(new_ind), n_unique_lbls(1));
+end
+
+new_lbl = valid_lbls{new_ind};
+% note: we remove any label fields below if they're not of the appropriate
+% length to correspond to the original parcellation
+kept_lbl_ind = find((n_unique_lbls <= n_unique_lbls(new_ind)) & ...
+    (n_lbls == length(obj.(labelfield))));
+
+rm_lbl_ind = find(~ismember(1:length(valid_lbls), 1:length(kept_lbl_ind)));
+
+% merge high resolution regions
+[new_lbl_parcels, lbl_exp] = unique(obj.(new_lbl),'stable');
+new_parcels = cell(1,length(new_lbl_parcels));
+
+% init progress watcher
+n_reg = length(new_lbl_parcels);
+last_msg = sprintf('Creating new region 1/%d',n_reg);
+fprintf('%s',last_msg);
+
+% working with sparse matrices is much faster for atlases with many
+% parcels, while atlases with few parcels are likely to run quickly
+% regardless
+obj.probability_maps = sparse(obj.probability_maps);
+for i = 1:length(new_lbl_parcels)
+    % this can be parallelized but it's very memory intensive for some
+    % atlases (e.g. canlab2023) and parallelization makes the problem worse
+    % by several factors.
+
+    % update progress watcher
+    delete_last_chars = length(last_msg);
+    fprintf('%s',char(8*ones(1,delete_last_chars)))
+    new_msg = sprintf('Creating new region %d/%d',i,n_reg);
+    fprintf('%s',new_msg);
+    last_msg = new_msg;
+
+    new_parcel{i} = obj.select_atlas_subset(new_lbl_parcels(i), new_lbl, 'exact', 'flatten');
+end
+fprintf('\n');
+
+% init progress watcher
+n_reg = length(new_parcel);
+last_msg = sprintf('Merging new region 1/%d',n_reg);
+fprintf('%s',last_msg);
+
+% Combine parcels
+new_atlas_obj = new_parcel{1};
+for i = 2:length(new_parcel)
+    % update progress watcher
+    delete_last_chars = length(last_msg);
+    fprintf('%s',char(8*ones(1,delete_last_chars)))
+    new_msg = sprintf('Merging new region %d/%d',i,n_reg);
+    fprintf('%s',new_msg);
+    last_msg = new_msg;
+
+    new_atlas_obj = new_atlas_obj.merge_atlases(new_parcel{i},'noverbose');
+end
+fprintf('\n');
+new_atlas_obj.probability_maps = sparse(new_atlas_obj.probability_maps);
+
+for i = 1:length(kept_lbl_ind)
+    % increment by 1 
+    old_lbls = obj.(valid_lbls{kept_lbl_ind(i)});
+
+    new_atlas_obj.(valid_lbls{i}) = old_lbls(lbl_exp(1:length(lbl_exp)));
+end
+
+% delete finer parcellations
+for ind = rm_lbl_ind(:)'
+    new_atlas_obj.(valid_lbls{ind}) = {};
+end
+
+new_atlas_obj.references = unique(new_atlas_obj.references,'rows');
+
+end

CanlabCore/@atlas/merge_atlases.m

@@ -40,6 +40,7 @@
 
             case 'always_replace', always_replace = true; varargin{i} = [];
             case 'noreplace', doreplacevoxels = false; varargin{i} = [];
+            case 'noverbose', doverbose = false; varargin{i} = [];
 
             otherwise , warning(['Unknown input string option:' varargin{i}]);
         end
@@ -50,11 +51,12 @@
 % Subset
 % -------------------------------------------------------------------------
 
-if ~isempty(varargin) && any(cellfun(@isvector, varargin) || cellfun(@iscell, varargin))
+if ~isempty(varargin) && any(cellfun(@isvector, varargin) | cellfun(@iscell, varargin))
 
     if doverbose, disp('Getting atlas subset with select_atlas_subset'); end
 
-    atlas_obj_to_add = select_atlas_subset(atlas_obj_to_add, varargin{:});
+    args = varargin(cellfun(@isvector, varargin) | cellfun(@iscell, varargin));
+    atlas_obj_to_add = select_atlas_subset(atlas_obj_to_add, args{:});
 
 end
 
@@ -135,7 +137,20 @@
         atlas_obj_to_add.probability_maps(wh, :) = 0;
     end
 
-    atlas_obj.probability_maps = [full(atlas_obj.probability_maps) full(atlas_obj_to_add.probability_maps)];
+    % concatenating sparse matrices first is much faster than converting to full and the
+    % concatenating full matrices so instead of doing a double conversion
+    % and concatenating, try concatenating sparse matricies if either of
+    % them are sparse and then converting to full
+    %atlas_obj.probability_maps = [full(atlas_obj.probability_maps) full(atlas_obj_to_add.probability_maps)];
+    if issparse(atlas_obj.probability_maps) || issparse(atlas_obj_to_add.probability_maps)
+        if ~issparse(atlas_obj.probability_maps)
+            atlas_obj.probability_maps = sparse(double(atlas_obj.probability_maps));
+        end
+        if ~issparse(atlas_obj_to_add.probability_maps)
+            atlas_obj_to_add.probability_maps = sparse(double(atlas_obj_to_add.probability_maps));
+        end
+    end
+    atlas_obj.probability_maps = full([atlas_obj.probability_maps atlas_obj_to_add.probability_maps]);
 
     atlas_obj = probability_maps_to_region_index(atlas_obj);
 

CanlabCore/@atlas/select_atlas_subset.m

@@ -14,8 +14,11 @@
 % Options:
 % 'flatten' : Flatten integer vector in .dat to a single 1/0 mask. Good for
 % creating masks from combinations of regions, or adding a set to another
-% atlas as a single atlas region. .probability_maps reduced to single max
-% map
+% atlas as a single atlas region. .probability_maps is reestimated as
+% posterior probability of union of constituent regions under conditional 
+% independence assumption. If probability maps don't sum to 1 across voxels
+% then conditional independence assumption is violated and we default to
+% using the maximum probability instead.
 %
 % 'labels_2' : If you enter any field name in the object, e.g., labels_2,
 % the function will search for keyword matches here instead of in obj.labels.
@@ -232,10 +235,25 @@
     % this is wrong. Consider, you could have a voxel that has 50% 
     % probability of being voxel A and 50% probability of being voxel B.
     % the max value is %50, but the probability that it's one or the other
-    % must be 100%. Let's try assuming conditional independence instead.
-    % There are a number of assumptions we need to make for a better
-    % solution though, so let this be the default
-    obj_subset.probability_maps = max(obj_subset.probability_maps, [], 2);
+    % must be 100%. Let's check if it's safe to assume conditional 
+    % independence instead.
+    %
+    % We don't know if the atlas labels are exhaustive, but if the sum of
+    % atlas label probabilities exceeds 1 we can conclude that P(A,B) ~=
+    % P(A) + P(B), because the total density must sum to 1 to be a valid
+    % probability measure. We check if probabilities exceed 1 with some
+    % tolerance to deal with floating point errors.
+    hasdata = any(obj.probability_maps,2);
+    tol = 1e-7;
+    if all(sum(obj.probability_maps(hasdata,:), 2) - 1 < tol)
+        obj_subset.probability_maps = sum(obj_subset.probability_maps, 2);
+    else
+        % we have evidence that probabilities are not conditionally
+        % independent
+        warning('Probabilities don''t sum to 1, so cannot assume condition independence. New map will use maximum probability instead of sum of probabilities.')
+        obj_subset.probability_maps = max(obj_subset.probability_maps, [], 2);
+    end
+
 
     % Add labels for combined mask
 

CanlabCore/@brainpathway/brainpathway.m

@@ -101,34 +101,21 @@
 % Creating class instances
 % -----------------------------------------------------------------------
 % You can create an empty object by using:
-% fmri_dat = fmri_data
-% - fmri_dat is the object.
-% - It will be created with a standard brain mask, brainmask.nii
-% - This image should be placed on your Matlab path
-% - The space information is stored in fmri_dat.volInfo
-% - Data is stored in fmri_dat.dat, in a [voxels x images] matrix
-% - You can replace or append data to the fmri_dat.dat field.
-%
-% You can create an object by assembling an image_vector object from parts
-% (entering fields) and converting using fmri_obj = fmri_data(image_vec_obj)
-%
-% You can create an fmri_data object with extacted image data.
-% - Let "imgs" be a string array or cell array of image names
-% - This command creates an object with your (4-D) image data:
-% - fmri_dat = fmri_data(imgs);
-% - Images can be zipped (.gz) as well. fmri_data() will unpack them.
-% - Only values in the standard brain mask, brainmask.nii, will be included.
-% - This saves memory by reducing the number of voxels saved dramatically.
-%
-% You can specify any mask you'd like to extract data from.
-% - Let "maskimagename" be a string array with a mask image name.
-% - this command creates the object with data saved in the mask:
-% - fmri_dat = fmri_data(imgs, maskimagename);
-% - The mask information is saved in fmri_dat.mask
-%
-% e.g., this extracts data from images within the standard brain mask:
-% dat = fmri_data(imgs, which('brainmask.nii'));
-%
+% b = brainpathways();
+% 
+% You can assign optional inputs with the following fields, followed by
+% data for the relevant field:
+% - any atlas_class object
+% 'voxel_dat'   
+% 'node_dat'    
+% 'region_dat'    
+% 'network_dat'    
+% 'partition_dat' 
+% 'node_weights'    
+% 'connectivity'    
+% 'connections_apriori'    
+% 'additional_info'
+%             
 % Defining the space of the extracted data
 % -----------------------------------------------------------------------
 % Note: There are two options for defining the space (i.e., coordinates/voxels)
@@ -285,7 +272,7 @@
             % allowable_inputs = properties(brainpathway)';       % should
             % be row cell vector of property names. This is recursive if we
             % do this though...
-            allowable_inputs = {'region_atlas'    'voxel_dat'    'node_dat'    'region_dat'    'network_dat'    'partition_dat' 'node_weights'    'connectivity'    'connections_apriori'    'additional_info'};
+            allowable_inputs = {'region_atlas'  'voxel_dat'    'node_dat'    'region_dat'    'network_dat'    'partition_dat' 'node_weights'    'connectivity'    'connections_apriori'    'additional_info'};
 
             % optional inputs with default values - each keyword entered will create a variable of the same name
 
@@ -327,7 +314,7 @@
 
                 switch class(varargin{i})
 
-                    case 'atlas', obj.region_atlas = varargin{i}; % input_atlas = true;
+                    case {'atlas'}, obj.region_atlas = varargin{i}; % input_atlas = true;
 
                     case 'char'
                         switch varargin{i}

CanlabCore/@fmri_data/plot.m

@@ -177,7 +177,7 @@
 
 allowable_inputs = {'plotmethod' 'doorthviews' 'dooutliers' 'dorunmeanmontages'};
 
-keyword_inputs = {'noorthviews' 'nooutliers' 'norunmontages' 'montage' 'montages' 'means_for_unique_Y'};
+keyword_inputs = {'noorthviews' 'nooutliers' 'norunmontages' 'nomontages' 'montage' 'montages' 'means_for_unique_Y'};
 
 % optional inputs with default values - each keyword entered will create a variable of the same name
 
@@ -210,7 +210,7 @@
 
             case 'nooutliers', dooutliers = false;
 
-            case 'norunmontages', dorunmeanmontages = false;
+            case {'norunmontages', 'nomontages'}, dorunmeanmontages = false;
 
             case {'montage' 'montages'}, plotmethod = 'montages'; % montages only
 
@@ -231,7 +231,7 @@
 switch plotmethod
     %  ==============================================================
     case 'data'
-        %  ==============================================================
+    %  ==============================================================
 
         if isempty(fmridat.dat)
             warning('No data in .dat field.');
@@ -385,7 +385,7 @@
         drawnow
 
         % ---------------------------------------------------------------
-        % Global mean vs. time
+        % Global mean vs. image series (i.e., time)
         % ---------------------------------------------------------------
         if size(fmridat.dat,2) > 1
 

CanlabCore/@fmri_data/rescale.m

@@ -28,6 +28,12 @@
 %                           *Note: these methods must exclude invalid (0 or
 %                           NaN) voxels image-wise. Some images (but not others) in an
 %                           object may be missing some voxels.
+%     - prctileimages       - Convert image values to percentile scores. This
+%                           normalizes the range and distribution of all
+%                           images. It is like ranking but MUCH faster for
+%                           images with many elements.
+%                           - Normalize each image by number of valid values, so images in set
+%                           will be on the same scale if some are missing voxels
 %
 %     - l2norm_images       divide each image by its l2 norm, multiply by sqrt(n valid voxels)
 %     - divide_by_csf_l2norm  divide each image by CSF l2 norm. requires MNI space images for ok results!
@@ -131,6 +137,47 @@
 
         fmridat.history{end+1} = 'Ranked images (columns) across voxels';
 
+
+    case 'prctileimages'
+
+        res = 0.1;  
+
+        newobj = fmridat;
+
+        for i = 1:size(fmridat.dat, 2)
+
+            d = fmridat.dat(:,i);
+
+            % Consider missing voxels, and exclude case-wise (image-wise)
+            ismissing = d == 0 | isnan(d);
+            d(ismissing) = NaN;
+
+            if ~all(d == 0 | isnan(d))
+
+                x = [-Inf prctile(d, 0:res:100)];
+
+                for j = 1:length(x) - 1
+
+                    wh = d <= x(j + 1) & d > x(j);
+                    newobj.dat(wh, i) = j;
+
+                end % bins
+
+                % normalize by number of valid values, so images in set
+                % will be on the same scale if some are missing voxels
+                d = d ./ length(d(~ismissing));
+
+            end % if image has values
+
+        end % images
+
+        newobj.dat(isnan(newobj.dat)) = 0; % Replace with 0 for compatibility with image format
+
+        fmridat = newobj;
+
+        fmridat.history{end+1} = 'Ranked images (columns) across voxels';
+
+
     case 'centerimages'
 
         % center images (observations)