From 1a729a9ee6ddb788f8b60421cb0fe8958e6b433e Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Wed, 2 Oct 2019 08:23:35 +0300 Subject: [PATCH 1/8] makeROIsExactlyContiguous.m: only goes through the whole loop if ROIs actually share voxels --- .../makeROIsExactlyContiguous.m | 32 ++++++++++--------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m b/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m index a455bdd8c..55886910f 100644 --- a/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m +++ b/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m @@ -36,23 +36,25 @@ coords1 = rois(iRoi).coords'; coords2 = rois(jRoi).coords'; [commonCoordinates, indexROI1, indexROI2] = intersect(coords1,coords2,'rows'); - %remove common coordinates from ROIs 1 and 2 - coords1 = setdiff(coords1,commonCoordinates,'rows'); - coords2 = setdiff(coords2,commonCoordinates,'rows'); - %attribute common coordinates to one or the other ROI depending on distance - belongsToROI1 = false(size(commonCoordinates,1),1); - for iCoords = 1:size(commonCoordinates,1) - %compute distance between these coordinates and all coordinates unique to either both ROI - distanceCoords1 = sqrt(sum((repmat(commonCoordinates(iCoords,1:3),size(coords1,1),1) - coords1(:,1:3)).^2,2)); - distanceCoords2 = sqrt(sum((repmat(commonCoordinates(iCoords,1:3),size(coords2,1),1) - coords2(:,1:3)).^2,2)); - %identify closest ROI - if min(distanceCoords1) < min(distanceCoords2) - belongsToROI1(iCoords) = true; + if ~isempty(commonCoordinates) + %remove common coordinates from ROIs 1 and 2 + coords1 = setdiff(coords1,commonCoordinates,'rows'); + coords2 = setdiff(coords2,commonCoordinates,'rows'); + %attribute common coordinates to one or the other ROI depending on distance + belongsToROI1 = false(size(commonCoordinates,1),1); + for iCoords = 1:size(commonCoordinates,1) + %compute distance between these coordinates and all coordinates unique to either both ROI + distanceCoords1 = sqrt(sum((repmat(commonCoordinates(iCoords,1:3),size(coords1,1),1) - coords1(:,1:3)).^2,2)); + distanceCoords2 = sqrt(sum((repmat(commonCoordinates(iCoords,1:3),size(coords2,1),1) - coords2(:,1:3)).^2,2)); + %identify closest ROI + if min(distanceCoords1) < min(distanceCoords2) + belongsToROI1(iCoords) = true; + end end + % delete coords that belong to the other ROI + rois(iRoi).coords(:,indexROI1(~belongsToROI1'))=[]; + rois(jRoi).coords(:,indexROI2(belongsToROI1'))=[]; end - % delete coords that belong to the other ROI - rois(iRoi).coords(:,indexROI1(~belongsToROI1'))=[]; - rois(jRoi).coords(:,indexROI2(belongsToROI1'))=[]; end end From bc58bfac643ac098d1d2dcc7f07d35ab755fa684 Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Wed, 2 Oct 2019 08:28:44 +0300 Subject: [PATCH 2/8] convertROICorticalDepth.m: removed mimimum possible value for minDepth (although mrParamsDialo does not seem to enforce this), moved removal of all voxels to the end of the function, and added some explanatory comments --- mrLoadRet/ROI/convertROICorticalDepth.m | 54 +++++++++++++------------ 1 file changed, 29 insertions(+), 25 deletions(-) diff --git a/mrLoadRet/ROI/convertROICorticalDepth.m b/mrLoadRet/ROI/convertROICorticalDepth.m index 632ebf4d8..8af4f682a 100644 --- a/mrLoadRet/ROI/convertROICorticalDepth.m +++ b/mrLoadRet/ROI/convertROICorticalDepth.m @@ -56,10 +56,10 @@ paramsInfo = {}; paramsInfo{end+1} = {'conversionType',{'Project through depth','Restrict to reference depth'},'type=popupmenu','If you set project through depth, then this will add all the voxels from each cortical depth that are in the same position as the ones at the reference depth. If you set to restrict to reference depth, this will remove any voxels that are not on the reference depth (note that you will still see some voxels on other depths, but those are voxels that exist at the reference depth--also, voxels that do not exist on this flat map will not be affected)'}; paramsInfo{end+1} = {'referenceDepth',referenceDepth,'min=0','max=1',incdecString,'The cortical depth to start from'}; - paramsInfo{end+1} = {'minDepth',minDepth,'min=0','max=1',incdecString,'The start depth'}; + paramsInfo{end+1} = {'minDepth',minDepth,'max=1',incdecString,'The minimum depth. Negative values will extend the ROI into white matter'}; paramsInfo{end+1} = {'depthStep',depthStep,'min=0','max=1',incdecString,'The depth step (i.e. we will go from minDepth:depthStep:maxDepth (skipping the reference depth), including or excluding voxels'}; - paramsInfo{end+1} = {'maxDepth',maxDepth,'min=0','max=1',incdecString,'The end depth'}; - paramsInfo{end+1} = {'excludeOtherVoxels',1,'type=checkbox','If ROI voxels exist oustide the projected surface, they will be remove. Uncheck to keep them. this option is ignored if restriction is selected'}; + paramsInfo{end+1} = {'maxDepth',maxDepth,'min=0','max=1',incdecString,'The maximum depth'}; + paramsInfo{end+1} = {'excludeOtherVoxels',1,'type=checkbox','If ROI voxels exist oustide the projected surface, they will be removed. Uncheck to keep them. This option is ignored if restriction is selected'}; if defaultParams params = mrParamsDefault(paramsInfo); else @@ -118,7 +118,7 @@ continue; end nVoxelsOriginalROI = size(roiBaseCoords,2); - % get base info + % get base info, including (rounded) base coordinates corresponding to the reference cortical depth baseVoxelSize = viewGet(v,'baseVoxelSize',baseNum); baseCoordMap = viewGet(v,'baseCoordMap',baseNum,params.referenceDepth); baseDims = baseCoordMap.dims; @@ -127,10 +127,12 @@ referenceBaseCoordMap = referenceBaseCoordMap(:); % get roi linear coordinates roiBaseCoordsLinear = mrSub2ind(baseDims,roiBaseCoords(1,:),roiBaseCoords(2,:),roiBaseCoords(3,:)); - % now find which baseCoords are in the current roi - [isInROI roiInBase] = ismember(referenceBaseCoordMap,roiBaseCoordsLinear); - % get the roi base coordinates that are found in base + % now find which baseCoords are in the current roi at the reference depth + [isInROI, roiInBase] = ismember(referenceBaseCoordMap,roiBaseCoordsLinear); + % get the roi base coordinates that are found in base at the reference depth roiInBase = unique(setdiff(roiInBase,0)); + % (note that here we could have used ismember(roiBaseCoordsLinear,referenceBaseCoordMap) instead, which is perhaps easier to understand) + % if we don't find most of the coordinates, then % probably good to complain and give up if (length(roiInBase)/length(roiBaseCoordsLinear)) < 0.1 @@ -140,49 +142,51 @@ end % make sure to keep the voxels at the reference depth roiBaseCoordsReferenceLinear = roiBaseCoordsLinear(ismember(roiBaseCoordsLinear,referenceBaseCoordMap)); + % (Note that we could have used roiBaseCoordsLinear(roiInBase) instead here) if strcmp(params.conversionType,'Project through depth') - %clear all voxels if we're not keeping voxels outside the projection - if params.excludeOtherVoxels - % remove everything from the ROI - roiBaseCoords(4,:) = 1; - v = modifyROI(v,roiBaseCoords,base2roi,baseVoxelSize,0); - end roiBaseCoordsLinear=[]; % now get each cortical depth, and add/remove voxels corticalDepths = params.minDepth:params.depthStep:params.maxDepth; + % (negative cortical depths mean that the flat/surface base (and subsequently the ROI) will be extended into white matter baseCoordMap = viewGet(v,'baseCoordMap',baseNum,corticalDepths); for iDepth = 1:size(baseCoordMap.coords,5) - % get the coordinates at this depth + % get the (rounded) base coordinates at this depth baseCoords = round(baseCoordMap.coords(:,:,:,:,iDepth)); baseCoords = mrSub2ind(baseDims,baseCoords(:,:,:,1),baseCoords(:,:,:,2),baseCoords(:,:,:,3)); baseCoords = baseCoords(:); - % add the coordinates to our list + % find the baseCoords that are in the ROI at this depth and add them to our list roiBaseCoordsLinear = union(roiBaseCoordsLinear,baseCoords(isInROI)); end roiBaseCoordsLinear = roiBaseCoordsLinear(~isnan(roiBaseCoordsLinear)); % now convert back to regular coords - roiBaseCoords = []; - [roiBaseCoords(1,:) roiBaseCoords(2,:) roiBaseCoords(3,:)] = ind2sub(baseDims,roiBaseCoordsLinear); - roiBaseCoords(4,:) = 1; - % add them to the ROI - v = modifyROI(v,roiBaseCoords,base2roi,baseVoxelSize,1); + additionalRoiBaseCoords = []; + [additionalRoiBaseCoords(1,:), additionalRoiBaseCoords(2,:), additionalRoiBaseCoords(3,:)] = ind2sub(baseDims,roiBaseCoordsLinear); + additionalRoiBaseCoords(4,:) = 1; + %clear all existing voxels if we're not keeping voxels outside the projection + if params.excludeOtherVoxels + % remove everything from the ROI + roiBaseCoords(4,:) = 1; + v = modifyROI(v,roiBaseCoords,base2roi,baseVoxelSize,0); + end + % add the projected voxels to the ROI + v = modifyROI(v,additionalRoiBaseCoords,base2roi,baseVoxelSize,1); else % get current coords curROICoords = viewGet(v,'roiCoords',roinum); % remove them from the ROI v = modifyROI(v,roiBaseCoords,base2roi,baseVoxelSize,0); % but make sure we have the voxels at the reference depth - roiBaseCoords = []; - [roiBaseCoords(1,:) roiBaseCoords(2,:) roiBaseCoords(3,:)] = ind2sub(baseDims,roiBaseCoordsReferenceLinear); - roiBaseCoords(4,:) = 1; - v = modifyROI(v,roiBaseCoords,base2roi,baseVoxelSize,1); + additionalRoiBaseCoords = []; + [additionalRoiBaseCoords(1,:), additionalRoiBaseCoords(2,:), additionalRoiBaseCoords(3,:)] = ind2sub(baseDims,roiBaseCoordsReferenceLinear); + additionalRoiBaseCoords(4,:) = 1; + v = modifyROI(v,additionalRoiBaseCoords,base2roi,baseVoxelSize,1); % and save for undo (note we do this instead of allowing % modifyROI to do it since we have called modifyROI twice) v = viewSet(v,'prevROIcoords',curROICoords); end disppercent(inf); - fprintf(1,'(convertROICorticalDepth) Number of voxels in original ROI: %d\t Number of voxels in modified ROI: %d\n',nVoxelsOriginalROI,size(roiBaseCoords,2)); + fprintf(1,'(convertROICorticalDepth) Number of voxels in original ROI: %d\t Number of voxels in modified ROI: %d\n',nVoxelsOriginalROI,size(additionalRoiBaseCoords,2)); end v = viewSet(v,'currentROI',currentROI); From 8aaf121d88f0155f66d4ee58cc723ef09a9e12da Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Wed, 2 Oct 2019 08:31:43 +0300 Subject: [PATCH 3/8] convertROICorticalDepth.m: started implementing option to exclude voxels that span two distant cortical locations (but not satisfactory yet) --- mrLoadRet/ROI/convertROICorticalDepth.m | 9 +++++++++ 1 file changed, 9 insertions(+) diff --git a/mrLoadRet/ROI/convertROICorticalDepth.m b/mrLoadRet/ROI/convertROICorticalDepth.m index 8af4f682a..f23ab86b4 100644 --- a/mrLoadRet/ROI/convertROICorticalDepth.m +++ b/mrLoadRet/ROI/convertROICorticalDepth.m @@ -60,6 +60,7 @@ paramsInfo{end+1} = {'depthStep',depthStep,'min=0','max=1',incdecString,'The depth step (i.e. we will go from minDepth:depthStep:maxDepth (skipping the reference depth), including or excluding voxels'}; paramsInfo{end+1} = {'maxDepth',maxDepth,'min=0','max=1',incdecString,'The maximum depth'}; paramsInfo{end+1} = {'excludeOtherVoxels',1,'type=checkbox','If ROI voxels exist oustide the projected surface, they will be removed. Uncheck to keep them. This option is ignored if restriction is selected'}; + paramsInfo{end+1} = {'allowProjectionThroughSulci',1,'type=checkbox','Voxels will be kept even if they also belong to another part of the cortical surface through a sulcus. Uncheck to exclude these voxels. Note that voxels projected to another part of the cortex through white matter (for instance in case minDepth is negative) will be excluded too. If the ROI is projected based on a flat map, only voxels on the flat map, not the whole surface, will be excluded. This option is ignored if restriction is selected'}; if defaultParams params = mrParamsDefault(paramsInfo); else @@ -159,6 +160,14 @@ roiBaseCoordsLinear = union(roiBaseCoordsLinear,baseCoords(isInROI)); end roiBaseCoordsLinear = roiBaseCoordsLinear(~isnan(roiBaseCoordsLinear)); + % exclude any projected voxel that ends up in a different part of cortex either through a sulcus or through white matter + if ~params.allowProjectionThroughSulci + % get the (rounded) base coordinates at all depths between 0 and 1 + baseCoords = round(baseCoordMap.coords(:,:,:,:,corticalDepths>=0 & corticalDepths<=1)); + baseCoords = mrSub2ind(baseDims,baseCoords(:,:,:,1,:),baseCoords(:,:,:,2,:),baseCoords(:,:,:,3,:)); + baseCoords = reshape(baseCoords,size(baseCoords,1)*size(baseCoords,2)*size(baseCoords,3),size(baseCoords,5)); + roiBaseCoordsLinear = setdiff(roiBaseCoordsLinear,baseCoords(~isInROI,:)); + end % now convert back to regular coords additionalRoiBaseCoords = []; [additionalRoiBaseCoords(1,:), additionalRoiBaseCoords(2,:), additionalRoiBaseCoords(3,:)] = ind2sub(baseDims,roiBaseCoordsLinear); From 497fbd5ae7fc47bb322dc23a8e89373420ec93ff Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Wed, 2 Oct 2019 14:51:40 +0300 Subject: [PATCH 4/8] projectCorticalDepth.m: finished implementing option to exclude ROI voxels belonging to distant locations on cortex --- mrLoadRet/ROI/convertROICorticalDepth.m | 98 +++++++++++++++++++------ 1 file changed, 77 insertions(+), 21 deletions(-) diff --git a/mrLoadRet/ROI/convertROICorticalDepth.m b/mrLoadRet/ROI/convertROICorticalDepth.m index f23ab86b4..7541197a6 100644 --- a/mrLoadRet/ROI/convertROICorticalDepth.m +++ b/mrLoadRet/ROI/convertROICorticalDepth.m @@ -32,6 +32,7 @@ if ieNotDefined('justGetParams'),justGetParams = 0;end if ieNotDefined('defaultParams'),defaultParams = 0;end +if ieNotDefined('distanceThreshold'), distanceThreshold = 2; end %distance threshold (in mm) to exclude % number of rois numrois = viewGet(v,'numberofrois'); @@ -41,6 +42,7 @@ end if ieNotDefined('params') + askForParams = 1; % get cortical depth corticalDepth = viewGet(v,'corticalDepth'); referenceDepth= mean(corticalDepth); @@ -53,26 +55,39 @@ end depthStep = 1/(mrGetPref('corticalDepthBins')-1); incdecString = sprintf('incdec=[-%f %f]',depthStep,depthStep); - paramsInfo = {}; - paramsInfo{end+1} = {'conversionType',{'Project through depth','Restrict to reference depth'},'type=popupmenu','If you set project through depth, then this will add all the voxels from each cortical depth that are in the same position as the ones at the reference depth. If you set to restrict to reference depth, this will remove any voxels that are not on the reference depth (note that you will still see some voxels on other depths, but those are voxels that exist at the reference depth--also, voxels that do not exist on this flat map will not be affected)'}; - paramsInfo{end+1} = {'referenceDepth',referenceDepth,'min=0','max=1',incdecString,'The cortical depth to start from'}; - paramsInfo{end+1} = {'minDepth',minDepth,'max=1',incdecString,'The minimum depth. Negative values will extend the ROI into white matter'}; - paramsInfo{end+1} = {'depthStep',depthStep,'min=0','max=1',incdecString,'The depth step (i.e. we will go from minDepth:depthStep:maxDepth (skipping the reference depth), including or excluding voxels'}; - paramsInfo{end+1} = {'maxDepth',maxDepth,'min=0','max=1',incdecString,'The maximum depth'}; - paramsInfo{end+1} = {'excludeOtherVoxels',1,'type=checkbox','If ROI voxels exist oustide the projected surface, they will be removed. Uncheck to keep them. This option is ignored if restriction is selected'}; - paramsInfo{end+1} = {'allowProjectionThroughSulci',1,'type=checkbox','Voxels will be kept even if they also belong to another part of the cortical surface through a sulcus. Uncheck to exclude these voxels. Note that voxels projected to another part of the cortex through white matter (for instance in case minDepth is negative) will be excluded too. If the ROI is projected based on a flat map, only voxels on the flat map, not the whole surface, will be excluded. This option is ignored if restriction is selected'}; - if defaultParams - params = mrParamsDefault(paramsInfo); - else - % put up some parameter choices - params = mrParamsDialog(paramsInfo,'ROI cortical depth conversion'); - end - % now select rois - % put up a dialog with rois to select - if defaultParams - params.roiList = viewGet(v,'curROI'); - else - params.roiList = selectInList(v,'rois'); + while askForParams + paramsInfo = {}; + paramsInfo{end+1} = {'conversionType',{'Project through depth','Restrict to reference depth'},'type=popupmenu','If you set project through depth, then this will add all the voxels from each cortical depth that are in the same position as the ones at the reference depth. If you set to restrict to reference depth, this will remove any voxels that are not on the reference depth (note that you will still see some voxels on other depths, but those are voxels that exist at the reference depth--also, voxels that do not exist on this flat map will not be affected)'}; + paramsInfo{end+1} = {'referenceDepth',referenceDepth,'min=0','max=1',incdecString,'The cortical depth to start from'}; + paramsInfo{end+1} = {'minDepth',minDepth,'max=1',incdecString,'The minimum depth. Negative values will extend the ROI into white matter'}; + paramsInfo{end+1} = {'depthStep',depthStep,'min=0','max=1',incdecString,'The depth step (i.e. we will go from minDepth:depthStep:maxDepth (skipping the reference depth), including or excluding voxels'}; + paramsInfo{end+1} = {'maxDepth',maxDepth,'min=0','max=1',incdecString,'The maximum depth'}; + paramsInfo{end+1} = {'excludeOtherVoxels',1,'type=checkbox','If ROI voxels exist oustide the projected surface, they will be removed. Uncheck to keep them. This option is ignored if restriction is selected'}; + paramsInfo{end+1} = {'allowProjectionThroughSulci',1,'type=checkbox','Voxels will be kept even if they also belong to another part of the cortical surface through a sulcus. Uncheck to exclude these voxels. Note that voxels projected to another part of the cortex through white matter (for instance in case minDepth is negative) will be excluded too. If the ROI is projected based on a flat map, only voxels on the flat map, not the whole surface, will be excluded. This option is ignored if restriction is selected'}; + if defaultParams + params = mrParamsDefault(paramsInfo); + else + % put up some parameter choices + params = mrParamsDialog(paramsInfo,'ROI cortical depth conversion'); + end + % Abort if params empty + if ieNotDefined('params'),return,end + + if 0 + %checks on params here if needed + else + askForParams = 0; + % now select rois + % put up a dialog with rois to select + if defaultParams + params.roiList = viewGet(v,'curROI'); + else + params.roiList = selectInList(v,'rois'); + if isempty(params.roiList) + askForParams = 1; + end + end + end end end if isempty(params),return,end @@ -109,6 +124,13 @@ baseNum = viewGet(v,'curBase'); end end + % check the base type to see if it's compatible with the current implementation of params.allowProjectionThroughSulci + if ~params.allowProjectionThroughSulci + if viewGet(v,'basetype',baseNum) == 2 + mrWarnDlg('(convertROICorticalDepth) Unchecking allowProjectionThroughSulci parameters is not yet implemented for surface bases') + return + end + end % get the roi transformation in order to set the coordinates later base2roi = viewGet(v,'base2roi',roinum,baseNum); % get the roiBaseCoords @@ -122,6 +144,7 @@ % get base info, including (rounded) base coordinates corresponding to the reference cortical depth baseVoxelSize = viewGet(v,'baseVoxelSize',baseNum); baseCoordMap = viewGet(v,'baseCoordMap',baseNum,params.referenceDepth); + mapDims = size(baseCoordMap.coords); baseDims = baseCoordMap.dims; baseCoordMap = round(baseCoordMap.coords); referenceBaseCoordMap = mrSub2ind(baseDims,baseCoordMap(:,:,:,1),baseCoordMap(:,:,:,2),baseCoordMap(:,:,:,3)); @@ -145,6 +168,38 @@ roiBaseCoordsReferenceLinear = roiBaseCoordsLinear(ismember(roiBaseCoordsLinear,referenceBaseCoordMap)); % (Note that we could have used roiBaseCoordsLinear(roiInBase) instead here) + % if excluding voxels that belong to two distant locations of the cortex, we need to compute the shortest distance of all elements + % of the flat map or surface to the ROI, in flat/surface space. (This is not yet implemented for surfaces and would require computing the Dijkstra distance) + if ~params.allowProjectionThroughSulci + [flatCoordsX, flatCoordsY] = meshgrid(1:mapDims(1),1:mapDims(2)); %compute the coordinates of the points on the flat map (i.e. in flat space). + % (Distance calculations could be done using the actual surface locations, but this would require computing Dijkstra distances. + % Easier like this and sufficient for our purposes, until the same is implemented for surfaces) + % find coordinates of the ROI on the flat map + flatCoordsRoiX = flatCoordsX(isInROI); + flatCoordsRoiY = flatCoordsY(isInROI); + % for each pixel of the flat map that corresponds to a base voxels, but that does not belong to the ROI, + % compute its shortest distance to any pixel in the ROI (in flat space) + minDistanceToROI = zeros(mapDims(1)*mapDims(2),1); + for iPixel = find(~isInROI & ~isnan(referenceBaseCoordMap))' + minDistanceToROI(iPixel) = min(sqrt((flatCoordsX(iPixel) - flatCoordsRoiX).^2 + (flatCoordsY(iPixel) - flatCoordsRoiY).^2)); + end + % now find flat map pixels that are a minimum distance from any pixel in the ROI. + % first compute approximate pixel size (based on surface coordinates only at the reference depth) + % separate x,y and z coordinates of flat map in base space and convert to mm + xBaseCoordMap = baseVoxelSize(1)*baseCoordMap(:,:,1,1); + yBaseCoordMap = baseVoxelSize(2)*baseCoordMap(:,:,1,2); + zBaseCoordMap = baseVoxelSize(3)*baseCoordMap(:,:,1,3); + % remove pixels that do not index a location on the surface + xBaseCoordMap(isnan(referenceBaseCoordMap))=NaN; + yBaseCoordMap(isnan(referenceBaseCoordMap))=NaN; + zBaseCoordMap(isnan(referenceBaseCoordMap))=NaN; + %compute pixel size in pixels + pixelSize(1) = nanmean(nanmean(sqrt(diff(xBaseCoordMap,1,1).^2 + diff(yBaseCoordMap,1,1).^2 + diff(zBaseCoordMap,1,1).^2))); + pixelSize(2) = nanmean(nanmean(sqrt(diff(xBaseCoordMap,1,2).^2 + diff(yBaseCoordMap,1,2).^2 + diff(zBaseCoordMap,1,2).^2))); + %compute distance threshold in pixels 2 based on approximate pixel size + isFarEnoughFromROI = minDistanceToROI > (distanceThreshold / min(pixelSize)); + end + if strcmp(params.conversionType,'Project through depth') roiBaseCoordsLinear=[]; % now get each cortical depth, and add/remove voxels @@ -166,7 +221,8 @@ baseCoords = round(baseCoordMap.coords(:,:,:,:,corticalDepths>=0 & corticalDepths<=1)); baseCoords = mrSub2ind(baseDims,baseCoords(:,:,:,1,:),baseCoords(:,:,:,2,:),baseCoords(:,:,:,3,:)); baseCoords = reshape(baseCoords,size(baseCoords,1)*size(baseCoords,2)*size(baseCoords,3),size(baseCoords,5)); - roiBaseCoordsLinear = setdiff(roiBaseCoordsLinear,baseCoords(~isInROI,:)); + % exclude any voxel that is also part of the flat map (or surface) at a location away form the ROI + roiBaseCoordsLinear = setdiff(roiBaseCoordsLinear,baseCoords(isFarEnoughFromROI,:)); end % now convert back to regular coords additionalRoiBaseCoords = []; From b2269de404ed7cba9c913503a321cea03ff96f11 Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Wed, 2 Oct 2019 23:02:32 +0300 Subject: [PATCH 5/8] makeROIsExactlyContiguous.m: handles ROIs that have non-rounded and non-unique voxels --- .../makeROIsExactlyContiguous.m | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) diff --git a/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m b/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m index 55886910f..e0384769a 100644 --- a/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m +++ b/mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m @@ -30,12 +30,12 @@ end end - for iRoi = 1:length(rois) for jRoi = iRoi+1:length(rois) - coords1 = rois(iRoi).coords'; - coords2 = rois(jRoi).coords'; - [commonCoordinates, indexROI1, indexROI2] = intersect(coords1,coords2,'rows'); + % first ensure that all voxel coordinates are rounded and unique + coords1 = unique(round(rois(iRoi).coords'),'rows'); + coords2 = unique(round(rois(jRoi).coords'),'rows'); + [commonCoordinates, indexROI1, indexROI2] = intersect(coords1,coords2,'rows'); % indexROI1, indexROI2 used to be used below if ~isempty(commonCoordinates) %remove common coordinates from ROIs 1 and 2 coords1 = setdiff(coords1,commonCoordinates,'rows'); @@ -51,9 +51,14 @@ belongsToROI1(iCoords) = true; end end - % delete coords that belong to the other ROI - rois(iRoi).coords(:,indexROI1(~belongsToROI1'))=[]; - rois(jRoi).coords(:,indexROI2(belongsToROI1'))=[]; +% % delete coords that belong to the other ROI +% rois(iRoi).coords(:,indexROI1(~belongsToROI1'))=[]; +% rois(jRoi).coords(:,indexROI2(belongsToROI1'))=[]; + % instead of deleting common voxels, replace all voxels in each ROI by its unique voxels + % and the common voxels that have been attributed to it + % (replacing is necessary because coordinates might have been rounded and duplciates removed) + rois(iRoi).coords = [coords1; commonCoordinates(belongsToROI1,:)]'; + rois(jRoi).coords = [coords2; commonCoordinates(~belongsToROI1,:)]'; end end end From 26931026a5579217c0b73d7f3237d3c76ec009bd Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Thu, 10 Oct 2019 20:24:56 +0300 Subject: [PATCH 6/8] fweAdjust.m: makes sure p values are in a column vector --- .../Plugin/GLM_v2/newGlmAnalysis/fweAdjust.m | 32 ++++++++++++++----- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/fweAdjust.m b/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/fweAdjust.m index 9b334e8c6..2556dec84 100644 --- a/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/fweAdjust.m +++ b/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/fweAdjust.m @@ -1,7 +1,7 @@ % fweAdjust.m % % $Id$ -% usage: adjustedP = fweAdjust(p,params,) +% usage: adjustedP = fweAdjust(p,params,) % by: julien besle, % date: 17/01/2011 % purpose: adjusts p-values for various familywise error control procedure @@ -14,6 +14,14 @@ function adjustedPdata = fweAdjust(p, params, numberTrueH0, lambda) +%this function assumes that p is a column vector +%if this is not the case, transpose +if size(p,1)==1 + p=p'; + transposed=true; +else + transposed=false; +end isNotNan = ~isnan(p); sizePdata = size(p); @@ -49,13 +57,17 @@ case 'Hommel' % Hommel (1988) Biometrika (1988), 75, 2, pp. 383-6 -%%% % p-adjustment algorithm provided by Wright (1992) -% % % % adjustedP=p; -% % % % for m=numberH0:-1:2 -% % % % cMin = min(m*p(numberH0-m+1:numberH0)./(m+(numberH0-m+1:numberH0)-numberH0)); -% % % % adjustedP(numberH0-m+1:numberH0) = max(adjustedP(numberH0-m+1:numberH0),cMin); -% % % % adjustedP(1:numberH0-m) = max(adjustedP(1:numberH0-m),min(cMin,m*p(1:numberH0-m))); -% % % % end + % p-adjustment algorithm provided by Wright (1992) +% % original version (note that this version assumes that p is a column vector) +% p=p'; +% adjustedP=p; +% for m=numberH0:-1:2 +% cMin = min(m*p(numberH0-m+1:numberH0)./(m+(numberH0-m+1:numberH0)-numberH0)); +% adjustedP(numberH0-m+1:numberH0) = max(adjustedP(numberH0-m+1:numberH0),cMin); +% adjustedP(1:numberH0-m) = max(adjustedP(1:numberH0-m),min(cMin,m*p(1:numberH0-m))); +% end +% p=p'; +% adjustedP = adjustedP'; %I think this version is clearer: adjustedP=p; @@ -73,3 +85,7 @@ adjustedP(sortingIndex) = adjustedP; adjustedPdata = NaN(sizePdata); adjustedPdata(isNotNan) = adjustedP; + +if transposed + adjustedPdata = adjustedPdata'; +end From 463f25f9cc1ca4036570eb38a7a070eacc9a91ae Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Mon, 14 Oct 2019 17:50:43 +0300 Subject: [PATCH 7/8] combineROIs.m: round coordinates to deal with ROIs that have non-integer coordinates --- mrLoadRet/ROI/combineROIs.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mrLoadRet/ROI/combineROIs.m b/mrLoadRet/ROI/combineROIs.m index b2721e3b8..858eeda79 100644 --- a/mrLoadRet/ROI/combineROIs.m +++ b/mrLoadRet/ROI/combineROIs.m @@ -73,7 +73,7 @@ elseif isempty(coords2) newCoords = coords2; else - newCoords = intersect(coords1,coords2,'rows'); + newCoords = intersect(round(coords1),round(coords2),'rows'); end case 'union' if isempty(coords1) @@ -89,7 +89,7 @@ elseif isempty(coords2) newCoords = coords1; else - newCoords = setxor(coords1,coords2,'rows'); + newCoords = setxor(round(coords1),round(coords2),'rows'); end case 'a not b' if isempty(coords1) @@ -97,7 +97,7 @@ elseif isempty(coords2) newCoords = coords1; else - newCoords = setdiff(coords1,coords2,'rows'); + newCoords = setdiff(round(coords1),round(coords2),'rows'); end otherwise error('unknown action: %s',action); From 469a91a606bc42cd548df77f338e96d25a5792b1 Mon Sep 17 00:00:00 2001 From: Julien Besle Date: Thu, 7 Nov 2019 03:52:58 +0200 Subject: [PATCH 8/8] (julien) added hrf model hrfCustom.m to add arbitrary HRF model using numerical values --- .../GLM_v2/newGlmAnalysis/glmAnalysisGUI.m | 2 +- .../Plugin/GLM_v2/newGlmAnalysis/hrfCustom.m | 63 +++++++++++++++++++ 2 files changed, 64 insertions(+), 1 deletion(-) create mode 100644 mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/hrfCustom.m diff --git a/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/glmAnalysisGUI.m b/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/glmAnalysisGUI.m index c32413da9..972992c87 100644 --- a/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/glmAnalysisGUI.m +++ b/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/glmAnalysisGUI.m @@ -94,7 +94,7 @@ askForParams = 1; % put group name on top of list to make it the default groupNames = putOnTopOfList(params.groupName,viewGet(thisView,'groupNames')); -hrfModelMenu = putOnTopOfList(params.hrfModel,{'hrfDoubleGamma','hrfFslFlobs','hrfDeconvolution','hrfBoxcar'}); +hrfModelMenu = putOnTopOfList(params.hrfModel,{'hrfDoubleGamma','hrfFslFlobs','hrfDeconvolution','hrfBoxcar','hrfCustom'}); analysisVolumeMenu = {'Whole volume'}; if nRois analysisVolumeMenu{end+1} = 'Loaded ROI(s)'; diff --git a/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/hrfCustom.m b/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/hrfCustom.m new file mode 100644 index 000000000..8f70a4c87 --- /dev/null +++ b/mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/hrfCustom.m @@ -0,0 +1,63 @@ +% hrfDeconvolution.m +% +% $Id$ +% usage: [params,hrf] = hrfCustom(params, sampleDuration, sampleDelay, defaultParams) +% by: julien besle +% date: 13/04/2010 +% purpose: returns the HRF specified as values in the parameters. If a time vector is specified +% checks that times correspond to actual TR and acquistion time (sampleDuration and sampleDelay) +% otherwise, assumes that HRF sample times correspond to those TR and acquisition times +% +function [params,hrf] = hrfCustom(params, sampleDuration,sampleDelay, defaultParams) + +if ~any(nargin == [1 2 3 4])% 5]) + help hrfCustom + return +end + +if ieNotDefined('defaultParams'),defaultParams = 0;end +if ieNotDefined('sampleDelay') + sampleDelay=sampleDuration/2; +end + +if ieNotDefined('params') + params = struct; +end +if fieldIsNotDefined(params,'description') + params.description = 'Custom HRF'; +end +if fieldIsNotDefined(params,'hrf') + [~, params.hrf] = hrfDoubleGamma([],sampleDuration,sampleDelay,1); + params.hrf = params.hrf'; +end +if fieldIsNotDefined(params,'hrfTimes') + params.hrfTimes = sampleDelay+sampleDuration*(0:length(params.hrf)-1); +end + +paramsInfo = {... + {'description', params.description, 'comment describing the hdr model'},... + {'hrf',params.hrf,'values of the the HRF'},... + {'hrfTimes',params.hrfTimes,'Times of the HRF samples'},... +}; + +if defaultParams + params = mrParamsDefault(paramsInfo); +else + params = mrParamsDialog(paramsInfo, 'Set Custom HRF parameters'); +end + +if nargout==1 + return +end + +%check that the times correspond to +if ~isequal(sampleDelay+sampleDuration*(0:length(params.hrf)-1), params.hrfTimes) + mrWarnDlg('(hrfCustom) HRF times are not compatible with TR and acquisition time'); + keyoard +else + if size(params.hrf,1)==1 + hrf = params.hrf'; + else + hrf = params.hrf; + end +end \ No newline at end of file