[ENH] refactor code for creation of ROIs from atlas

demos/atlas/create_roi_from_atlas.m

@@ -0,0 +1,19 @@
+opt.roi.atlas = 'wang';
+opt.roi.name = {'V1v', 'V1d'};
+opt.roi.dir = fullfile(pwd, 'derivatives', 'cpp_roi', 'group');
+
+spm_mkdir(opt.roi.dir);
+
+hemi = {'lh', 'rh'};
+
+for iHemi = 1:numel(hemi)
+
+  for iROI = 1:numel(opt.roi.name)
+
+    roiName = opt.roi.name{iROI};
+
+    imageName = extractRoiFromAtlas(opt.roi.dir, opt.roi.atlas, roiName, hemi{iHemi});
+
+  end
+
+end

demos/roi/label_and_extract_clusters.m

@@ -6,6 +6,8 @@
 % where each cluster has one label.
 %
 
+run ../../initCppRoi;
+
 gunzip(fullfile('inputs', '*.gz'));
 
 zMap = fullfile(pwd, 'inputs', 'visual motion_association-test_z_FDR_0.01.nii');

demos/roi/other_demo.m

@@ -1,5 +1,7 @@
 % (C) Copyright 2021 CPP ROI developers
 
+run ../../initCppRoi;
+
 gunzip(fullfile('inputs', '*.gz'));
 zMap = fullfile(pwd, 'inputs', 'visual motion_association-test_z_FDR_0.01.nii');
 

demos/roi/roi_script.m

@@ -10,6 +10,8 @@
 clear;
 clc;
 
+run ../../initCppRoi;
+
 %% ASSUMPTION
 %
 % This assumes that the 2 immages are in the same space (MNI, individual)

src/atlas/extractRoiFromAtlas.m

@@ -1,7 +1,37 @@
 % (C) Copyright 2021 CPP ROI developers
 
-function [imageName, imageVolume] = extractRoiFromAtlas(atlas, label, hemisphere)
+function roiImage = extractRoiFromAtlas(roiDir, atlas, roiName, hemisphere)
 
-  [imageName, imageVolume] = extractRoiFromLabel(sourceImage, label, hemisphere);
+  if strcmp(atlas, 'wang')
 
+    [maxProbaFiles, roiLabels] = getRetinoProbaAtlas();
+
+    if strcmp(hemisphere, 'lh')
+      sourceImage = maxProbaFiles(1, :);
+    else
+      sourceImage = maxProbaFiles(2, :);
+    end
+
+  end
+
+  roiIdx = strcmp(roiName, roiLabels.ROI);
+  label = roiLabels.label(roiIdx);
+
+  labelStruct = struct('ROI', [hemisphere roiName], ...
+                       'label', label);
+
+  roiImage = extractRoiByLabel(sourceImage, labelStruct);
+
+  nameStructure = struct( ...
+                         'space', 'MNI', ...
+                         'hemi', hemisphere, ...
+                         'desc', atlas, ...
+                         'label', roiName, ...
+                         'type', 'mask', ...
+                         'ext', '.nii');
+  newName = createFilename(nameStructure);
+
+  movefile(roiImage, fullfile(roiDir, newName));
+
+  roiImage = fullfile(roiDir, newName);
 end

src/atlas/getRetinoProbaAtlas.m

@@ -15,54 +15,13 @@
   %   PMID: 25452571
   %   Probabilistic Maps of Visual Topography in Human Cortex
   %
-  % If the data is not present in the ``cpp_spm/atlas/ProbAtlas_v4`` of the repo, it will be
-  % downloaded and unzipped
 
-  README_URL = 'http://scholar.princeton.edu/sites/default/files/napl/files/readme.txt';
-  ATLAS_URL = 'http://scholar.princeton.edu/sites/default/files/napl/files/probatlas_v4.zip';
+  unzipAtlas('wang');
 
-  ATLAS_FOLDER = fullfile( ...
-                          fileparts(mfilename('fullpath')), ...
-                          '..', '..',  'atlas');
+  atlasDir = returnAtlasDir('wang');
 
-  if ~exist(fullfile(ATLAS_FOLDER, 'ProbAtlas_v4'), 'dir')
+  maxProbaFiles = spm_select('FPList', fullfile(atlasDir, 'subj_vol_all'), '^.*_dseg.nii$');
 
-    mkdir(ATLAS_FOLDER);
-
-    try
-      urlwrite(ATLAS_URL, 'probatlas_v4.zip');
-      unzip('probatlas_v4.zip', ATLAS_FOLDER);
-    catch
-      system(sprintf('curl -L %s -o probatlas_v4.zip', ATLAS_URL));
-      unzip('probatlas_v4.zip', ATLAS_FOLDER);
-    end
-
-  end
-
-  %     urlwrite(README_URL, fullfile(ATLAS_FOLDER, 'ProbAtlas_v4', 'README.txt'));
-  %     urlread(README_URL)
-
-  maxProbaFiles = spm_select('FPListRec', ...
-                             fullfile(ATLAS_FOLDER, 'ProbAtlas_v4'), ...
-                             '^max.*.nii$');
-
-  if size(maxProbaFiles, 1) < 2
-
-    gunzip(fullfile(pwd, 'atlas', 'ProbAtlas_v4', 'subj_vol_all', 'max*.nii.gz'));
-
-    maxProbaFiles = spm_select('FPListRec', ...
-                               fullfile(ATLAS_FOLDER, 'ProbAtlas_v4'), ...
-                               '^max.*.nii$');
-
-  end
-
-  if size(maxProbaFiles, 1) < 2
-    error('no atlas present');
-  end
-
-  % the label file was created manually to be easier to load into SPM
-  roiLabels = spm_load(spm_select('FPListRec', ...
-                                  fullfile(ATLAS_FOLDER), ...
-                                  '^ROI_labels_ProbAtlas_v4.csv$'));
+  roiLabels = getRoiLabelLookUpTable('wang');
 
 end

src/atlas/getRoiLabelLookUpTable.m

@@ -7,26 +7,19 @@
     return
   end
 
+  atlasDir = returnAtlasDir(atlas);
+
   switch lower(atlas)
 
     case 'wang'
 
-      atlasDir = fullfile(fileparts(mfilename('fullpath')), ...
-                          '..', ...
-                          '..', ...
-                          'atlas');
-
-      roiLabelLUT = spm_load(fullfile(atlasDir, ...
-                                      'visual_topography_probability_atlas', ...
-                                      'LUT.csv'));
+      roiLabelLUT = spm_load(fullfile(atlasDir, 'LUT.csv'));
 
     case 'anatomy_toobox'
 
       anat_tb_URL = 'https://www.fz-juelich.de/SharedDocs/Downloads/INM/INM-1/DE/Toolbox/Toolbox_22c.html';
 
-      spmDir = spm('dir');
-
-      roiLabelLUT = fullfile(spmDir, 'toolbox', 'Anatomy', 'Anatomy_v22c_MPM.txt');
+      roiLabelLUT = fullfile(atlasDir, 'Anatomy_v22c_MPM.txt');
 
       if ~exist(roiLabelLUT, 'file')
         error('Did you install the spm Anatomy toolbox?\n\nDownload it from: %s', ...

src/atlas/returnAtlasDir.m

@@ -0,0 +1,23 @@
+% (C) Copyright 2021 CPP ROI developers
+
+function atlasDir = returnAtlasDir(atlas)
+
+  atlasDir = fullfile(fileparts(mfilename('fullpath')), '..', '..', 'atlas');
+
+  if nargin > 0
+
+    switch atlas
+
+      case 'wang'
+        atlasDir = fullfile(atlasDir, 'visual_topography_probability_atlas');
+
+      case 'anatomy_toobox'
+        atlasDir = fullfile(spm('dir'), 'toolbox', 'Anatomy');
+
+    end
+
+  end
+
+  atlasDir = spm_file(atlasDir, 'cpath');
+
+end

src/atlas/unzipAtlas.m

@@ -1,19 +1,27 @@
 % (C) Copyright 2021 CPP ROI developers
 
-function unzipAtlas()
+function unzipAtlas(atlas)
 
-  atlasDir = fullfile(fileparts(mfilename('fullpath')), '..', '..', 'atlas');
+  atlasDir = returnAtlasDir();
 
-  file = fullfile(atlasDir, 'visual_topography_probability_atlas.zip');
+  if strcmp(atlas, 'wang')
 
-  unzip(file, fileparts(file));
+    file = fullfile(atlasDir, 'visual_topography_probability_atlas.zip');
 
-  labelImages =  spm_select('FPList', ...
-                            fullfile(atlasDir, ...
-                                     'visual_topography_probability_atlas', ...
-                                     'subj_vol_all'), ...
-                            '^maxprob_vol_.*h.nii.gz$');
+    if ~exist(fullfile(atlasDir, 'visual_topography_probability_atlas'), 'dir')
 
-  gunzip(cellstr(labelImages));
+      unzip(file, fileparts(file));
+
+      labelImages =  spm_select('FPList', ...
+                                fullfile(atlasDir, ...
+                                         'visual_topography_probability_atlas', ...
+                                         'subj_vol_all'), ...
+                                '^.*_dseg.nii.gz$');
+
+      gunzip(cellstr(labelImages));
+
+    end
+
+  end
 
 end

src/roi/createRoi.m

@@ -158,12 +158,12 @@
         roiImage = specification.mask1;
         sphere = specification.mask2;
       end
-      
+
       % check that input image has at least enough voxels to include
       maskVol = spm_read_vols(spm_vol(roiImage));
       totalNbVoxels = sum(maskVol(:));
       if sphere.maxNbVoxels > totalNbVoxels
-          error('Number of voxels requested greater than the total number of voxels in this mask');
+        error('Number of voxels requested greater than the total number of voxels in this mask');
       end
 
       specification  = struct( ...
@@ -174,33 +174,33 @@
       hdr = spm_vol(roiImage);
       dim = diag(hdr.mat);
       radiusStep = min(abs(dim(1:3)));
-      
+
       % determine maximum radius to expand to
       maxRadius = hdr.dim .* dim(1:3)';
       maxRadius = max(abs(maxRadius));
 
-      fprintf(1, '\n Expansion:')
-      
+      fprintf(1, '\n Expansion:');
+
       while  true
         mask = createRoi('intersection', specification);
         mask.roi.radius = specification.mask2.radius;
-        
+
         fprintf(1, '\n radius: %0.2f mm; roi size: %i voxels', ...
-            mask.roi.radius, ...
-            mask.roi.size)
-        
+                mask.roi.radius, ...
+                mask.roi.size);
+
         if mask.roi.size > sphere.maxNbVoxels
           break
         end
-        
+
         if mask.roi.radius > maxRadius
-            error('sphere expanded beyond the dimension of the mask.')
+          error('sphere expanded beyond the dimension of the mask.');
         end
-        
+
         specification.mask2.radius = specification.mask2.radius + radiusStep;
       end
-      
-      fprintf(1, '\n')
+
+      fprintf(1, '\n');
 
       mask.xyz = sphere.location;
 
@@ -310,12 +310,12 @@
           fullfile(outputDir, tempFile));
 
   outputFile = fullfile(outputDir, roiName);
-      
+
   mars_rois2img(fullfile(outputDir, tempFile), ...
                 outputFile, ...
                 spm_vol(volumeDefiningImage));
   delete(fullfile(outputDir, tempFile));
-  
+
   % delete label files
   delete(fullfile(outputDir, '*_mask_labels.mat'));
 

tests/test_unit_convertToValidCamelCase.m

@@ -1,4 +1,4 @@
-% (C) Copyright 2021 CPP BIDS SPM-pipeline developers
+% (C) Copyright 2021 CPP ROI developers
 
 function test_suite = test_unit_convertToValidCamelCase %#ok<*STOUT>
   try % assignment of 'localfunctions' is necessary in Matlab >= 2016