Merge pull request #1 from julienbesle/current

updating fork

mrLoadRet/Plugin/GLM_v2/newGlmAnalysis/fweAdjust.m

@@ -1,7 +1,7 @@
 % fweAdjust.m
 %
 % $Id$
-% usage: adjustedP = fweAdjust(p,params,)
+% usage: adjustedP = fweAdjust(p,params,<numberTrueH0, lambda>)
 % 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

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)';

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

mrLoadRet/Plugin/GLM_v2/transformROIFunctions/makeROIsExactlyContiguous.m

@@ -30,29 +30,36 @@
  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');
- %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;
+ % 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');
+ 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'))=[];
+ % 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
- % delete coords that belong to the other ROI
- rois(iRoi).coords(:,indexROI1(~belongsToROI1'))=[];
- rois(jRoi).coords(:,indexROI2(belongsToROI1'))=[];
  end
 end
 

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,15 +89,15 @@
  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)
  newCoords = coords1;
  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);