Add fmri_surface_data object for CIFTI/GIFTI surface & grayordinate data

CanlabCore/@fmri_surface_data/apply_mask.m

@@ -0,0 +1,55 @@
+function obj = apply_mask(obj, mask, varargin)
+% apply_mask Mask a grayordinate object by zeroing out-of-mask grayordinates.
+%
+% :Usage:
+% ::
+%     obj = apply_mask(obj, mask)
+%     obj = apply_mask(obj, mask, 'invert')
+%
+% Surface analogue of fmri_data.apply_mask, greatly simplified per design
+% decision D5b: because all objects share a standardized grayordinate space,
+% masking is intrinsic -- there is no fmri_mask_image, no resampling, and no
+% empty-removal. Out-of-mask grayordinate rows are simply set to 0; .dat keeps
+% its full size and the geometry is unchanged.
+%
+% :Inputs:
+%   **obj:**  an fmri_surface_data object.
+%   **mask:** one of
+%       - logical/numeric vector [nGrayordinates x 1] (nonzero = keep)
+%       - another fmri_surface_data on the same space (nonzero/non-NaN = keep)
+%
+% :Optional Inputs:
+%   **'invert':** keep the complement (zero the in-mask grayordinates instead).
+%
+% :Outputs:
+%   **obj:** masked object (out-of-mask rows zeroed).
+%
+% :See also: fmri_surface_data, compare_space
+
+invert = any(strcmpi(varargin, 'invert'));
+
+if isa(mask, 'fmri_surface_data')
+    if compare_space(obj, mask) ~= 0
+        error('fmri_surface_data:apply_mask:space', ...
+            ['Mask is not on the same grayordinate space as the data ' ...
+             '(compare_space ~= 0). Resample first.']);
+    end
+    keep = any(mask.dat ~= 0 & ~isnan(mask.dat), 2);
+elseif islogical(mask) || isnumeric(mask)
+    keep = logical(mask(:));
+    if numel(keep) ~= size(obj.dat, 1)
+        error('fmri_surface_data:apply_mask:length', ...
+            'Mask vector length (%d) must equal the number of grayordinates (%d).', ...
+            numel(keep), size(obj.dat, 1));
+    end
+else
+    error('fmri_surface_data:apply_mask:type', ...
+        'mask must be a logical/numeric vector or an fmri_surface_data object.');
+end
+
+if invert, keep = ~keep; end
+
+obj.dat(~keep, :) = 0;
+obj.history{end+1} = sprintf('apply_mask: kept %d / %d grayordinates (zeroed the rest)', ...
+    nnz(keep), numel(keep));
+end

CanlabCore/@fmri_surface_data/apply_parcellation.m

@@ -0,0 +1,138 @@
+function [parcel_means, parcel_labels, parcel_table] = apply_parcellation(obj, parcels, varargin)
+% apply_parcellation Average grayordinate data within parcels of a surface atlas.
+%
+% :Usage:
+% ::
+%     parcel_means = apply_parcellation(obj, parcels)
+%     [parcel_means, labels, tbl] = apply_parcellation(obj, parcels, 'area')
+%
+% Computes the mean of each map over the grayordinates in each parcel of a
+% surface/grayordinate parcellation, mirroring image_vector.apply_parcellation.
+% Parcels with key 0 (and NaN) are treated as background / medial wall and are
+% excluded. Operates directly on .dat (no resampling needed when on the same
+% space).
+%
+% :Inputs:
+%   **obj:**     an fmri_surface_data object ([nGrayordinates x nMaps]).
+%   **parcels:** the parcellation, one of
+%       - an fmri_surface_data with integer label keys (e.g. a `.dlabel`),
+%         on the SAME grayordinate space (compare_space == 0); its label_table
+%         provides parcel names.
+%       - an integer vector [nGrayordinates x 1] of label keys.
+%
+% :Optional Inputs:
+%   **'area':** area-weight the average using per-vertex surface area (cortex)
+%               instead of an unweighted mean. (Subcortical voxels use unit weight.)
+%
+% :Outputs:
+%   **parcel_means:** [nMaps x nParcels] mean value per parcel.
+%   **parcel_labels:** 1 x nParcels cell of parcel names (or 'parcel_<key>').
+%   **parcel_table:**  table with columns key, label, n_grayordinates, (area).
+%
+% :Examples:
+% ::
+%     atl = fmri_surface_data(which('Gordon333.32k_fs_LR_Tian_Subcortex_S2.dlabel.nii'));
+%     s   = fmri_surface_data(which('transcriptomic_gradients.dscalar.nii'));
+%     pm  = apply_parcellation(s, atl);          % [nMaps x nParcels]
+%
+% :See also: fmri_surface_data, reparse_contiguous, surface_region, condf2indic
+
+use_area = any(strcmpi(varargin, 'area'));
+
+% ---- Resolve parcel keys + names ----
+names_by_key = containers.Map('KeyType', 'double', 'ValueType', 'char');
+if isa(parcels, 'fmri_surface_data')
+    if compare_space(obj, parcels) ~= 0
+        error('fmri_surface_data:apply_parcellation:space', ...
+            'Parcellation is not on the same grayordinate space (compare_space ~= 0).');
+    end
+    keys = round(double(parcels.dat(:, 1)));
+    if ~isempty(parcels.label_table)
+        for i = 1:numel(parcels.label_table)
+            names_by_key(parcels.label_table(i).key) = parcels.label_table(i).name;
+        end
+    end
+elseif isnumeric(parcels)
+    keys = round(double(parcels(:)));
+    if numel(keys) ~= size(obj.dat, 1)
+        error('fmri_surface_data:apply_parcellation:length', ...
+            'Parcel key vector length (%d) must equal the number of grayordinates (%d).', ...
+            numel(keys), size(obj.dat,1));
+    end
+else
+    error('fmri_surface_data:apply_parcellation:type', ...
+        'parcels must be an fmri_surface_data or an integer vector.');
+end
+
+% ---- Indicator matrix over positive keys (exclude 0 / NaN) ----
+ukeys = unique(keys(keys > 0 & ~isnan(keys)));
+nP = numel(ukeys);
+if nP == 0
+    error('fmri_surface_data:apply_parcellation:noparcels', 'No positive parcel keys found.');
+end
+indic = double(keys == ukeys');                 % [nGray x nP] 0/1
+
+% ---- Per-grayordinate weights ----
+if use_area
+    w = local_vertex_areas(obj);                % [nGray x 1]
+else
+    w = ones(size(obj.dat, 1), 1);
+end
+
+D = double(obj.dat);                            % [nGray x nMaps]
+wsum = (w' * indic);                            % [1 x nP] total weight per parcel
+sums = (D .* w)' * indic;                       % [nMaps x nP]
+parcel_means = sums ./ wsum;                    % [nMaps x nP]
+
+% ---- Labels + table ----
+parcel_labels = cell(1, nP);
+counts = sum(indic, 1);                         % grayordinates per parcel
+for i = 1:nP
+    if names_by_key.isKey(ukeys(i)) && ~isempty(names_by_key(ukeys(i)))
+        parcel_labels{i} = names_by_key(ukeys(i));
+    else
+        parcel_labels{i} = sprintf('parcel_%d', ukeys(i));
+    end
+end
+
+if nargout >= 3
+    parcel_table = table(ukeys(:), parcel_labels(:), counts(:), wsum(:), ...
+        'VariableNames', {'key', 'label', 'n_grayordinates', 'total_weight'});
+end
+end
+
+
+% -------------------------------------------------------------------------
+function w = local_vertex_areas(obj)
+% Per-grayordinate weight = barycentric vertex area on the cortical mesh
+% (1/3 of incident face areas); subcortical voxels get unit weight.
+w = ones(size(obj.dat, 1), 1);
+for mi = 1:numel(obj.brain_model.models)
+    m = obj.brain_model.models{mi};
+    rows = (m.start:(m.start + m.count - 1))';
+    if ~strcmp(m.type, 'surf'), continue; end
+    [V, F] = local_hemi_VF(obj.surface_space, m.struct);
+    va = local_face_vertex_area(V, F);          % [numvert x 1]
+    w(rows) = va(m.vertlist + 1);
+end
+end
+
+
+function [V, F] = local_hemi_VF(space, structname)
+% Prefer midthickness (true areas); fall back to inflated if unavailable.
+try
+    g = load_surface_geom(space, 'midthickness');
+catch
+    g = load_surface_geom(space, 'inflated');
+end
+if contains(upper(structname), 'LEFT'), V = g.vertices_lh; F = g.faces_lh;
+else, V = g.vertices_rh; F = g.faces_rh; end
+end
+
+
+function va = local_face_vertex_area(V, F)
+% Barycentric per-vertex area: each vertex gets 1/3 of each incident face's area.
+v1 = V(F(:,1), :); v2 = V(F(:,2), :); v3 = V(F(:,3), :);
+fa = 0.5 * sqrt(sum(cross(v2 - v1, v3 - v1, 2).^2, 2));   % [nFaces x 1]
+va = accumarray(F(:), repmat(fa, 3, 1) / 3, [size(V,1) 1]);
+end

CanlabCore/@fmri_surface_data/cat.m

@@ -0,0 +1,74 @@
+function obj = cat(obj, varargin)
+% cat Concatenate fmri_surface_data objects along the image (map) dimension.
+%
+% :Usage:
+% ::
+%     combined = cat(obj1, obj2, obj3, ...)
+%
+% Horizontally concatenates the [grayordinates x maps] data of two or more
+% fmri_surface_data objects (e.g. combining subjects/contrasts into one dataset),
+% along with their per-map annotations (X, Y, covariates, image_names,
+% metadata_table). All objects must share the same grayordinate space
+% (compare_space == 0); otherwise an error is raised. Surface analogue of
+% fmri_data.cat. No replace_empty step is needed (.dat is always full, D5b).
+%
+% :Inputs:
+%   **obj, varargin:** two or more fmri_surface_data objects on the same space.
+%
+% :Outputs:
+%   **obj:** a single fmri_surface_data with all maps concatenated.
+%
+% :Examples:
+% ::
+%     all_subs = cat(sub1, sub2, sub3);
+%     t = ttest(all_subs);
+%
+% :See also: horzcat, compare_space, fmri_surface_data
+
+for i = 1:numel(varargin)
+    o2 = varargin{i};
+    if ~isa(o2, 'fmri_surface_data')
+        error('fmri_surface_data:cat:type', 'All inputs must be fmri_surface_data objects.');
+    end
+    if compare_space(obj, o2) ~= 0
+        error('fmri_surface_data:cat:space', ...
+            ['Objects are not on the same grayordinate space (compare_space ~= 0). ' ...
+             'Resample to a common space before concatenating.']);
+    end
+
+    obj.dat = [obj.dat, o2.dat];
+
+    obj.X                = local_vcat(obj.X, o2.X);
+    obj.Y                = local_vcat(obj.Y, o2.Y);
+    obj.covariates       = local_vcat(obj.covariates, o2.covariates);
+    obj.image_names      = local_catcell(obj.image_names, o2.image_names);
+    obj.metadata_table   = local_cattable(obj.metadata_table, o2.metadata_table);
+    if ~isempty(obj.images_per_session) || ~isempty(o2.images_per_session)
+        obj.images_per_session = [obj.images_per_session(:); o2.images_per_session(:)]';
+    end
+end
+
+obj.removed_images = false(size(obj.dat, 2), 1);
+obj.removed_voxels = false(size(obj.dat, 1), 1);
+obj.history{end+1} = sprintf('cat: concatenated to %d maps', size(obj.dat, 2));
+end
+
+
+% -------------------------------------------------------------------------
+function v = local_vcat(a, b)
+if isempty(a) && isempty(b), v = []; return; end
+v = [a; b];          % per-map (rows = maps) fields
+end
+
+function c = local_catcell(a, b)
+if isempty(a) && isempty(b), c = {}; return; end
+a = reshape(a, [], 1); b = reshape(b, [], 1);
+c = [a; b];
+end
+
+function t = local_cattable(a, b)
+if isempty(a) && isempty(b), t = []; return; end
+if isempty(a), t = b; return; end
+if isempty(b), t = a; return; end
+t = [a; b];
+end

CanlabCore/@fmri_surface_data/compare_space.m

@@ -0,0 +1,66 @@
+function isdiff = compare_space(obj, obj2)
+% compare_space Compare the grayordinate spaces of two fmri_surface_data objects.
+%
+% :Usage:
+% ::
+%     isdiff = compare_space(obj, obj2)
+%
+% Surface analogue of image_vector.compare_space, preserving the same integer
+% return contract so callers (cat, apply_mask) can branch on it. Compares the
+% surface_space tag and the brain_model layout (per-model structure, count, and
+% surface vertex count), then the in-data selection (vertlist / voxlist).
+%
+% :Outputs:
+%   **isdiff:** integer code:
+%       - 0  same space and same in-data grayordinates
+%       - 1  different spaces (tag or model layout differ)
+%       - 2  no brain_model for one or more objects
+%       - 3  same space, but different in-data grayordinates (vertlist/voxlist
+%            or row count differ)
+%
+% :See also: image_vector.compare_space, resample_space, cat, apply_mask
+
+if isempty(obj.brain_model) || isempty(obj2.brain_model)
+    isdiff = 2;
+    return
+end
+
+% Space tag
+if ~strcmp(char(obj.surface_space), char(obj2.surface_space))
+    isdiff = 1;
+    return
+end
+
+m1 = obj.brain_model.models;
+m2 = obj2.brain_model.models;
+if numel(m1) ~= numel(m2)
+    isdiff = 1;
+    return
+end
+
+% Model layout: structure name, in-data count, and surface vertex count
+for i = 1:numel(m1)
+    same = strcmp(m1{i}.struct, m2{i}.struct) ...
+        && isequaln(m1{i}.numvert, m2{i}.numvert) ...
+        && strcmp(m1{i}.type, m2{i}.type);
+    if ~same
+        isdiff = 1;
+        return
+    end
+end
+
+% Same layout: now check the exact in-data selection
+isdiff = 0;
+if size(obj.dat, 1) ~= size(obj2.dat, 1)
+    isdiff = 3;
+    return
+end
+for i = 1:numel(m1)
+    if m1{i}.count ~= m2{i}.count ...
+            || ~isequal(m1{i}.vertlist, m2{i}.vertlist) ...
+            || ~isequal(m1{i}.voxlist, m2{i}.voxlist)
+        isdiff = 3;
+        return
+    end
+end
+end