Commit db4e8a55 authored by Gholamreza Salimi Khorshidi's avatar Gholamreza Salimi Khorshidi Committed by Duncan Mortimer
Browse files

Import of 1.05c sources

parents
*.nii.gz filter=lfs diff=lfs merge=lfs -text
*.RData filter=lfs diff=lfs merge=lfs -text
-------------------------------------------------------
FIX - FMRIB's ICA-based Xnoiseifier
Gholamreza Salimi-Khorshidi and Stephen Smith, FMRIB Analysis Group
Copyright (C) 2012-2013 University of Oxford
-------------------------------------------------------
SETUP
Requires: R, with the following commands run (maybe as root) to install packages:
install.packages("kernlab")
install.packages("ROCR")
install.packages("class")
install.packages("party")
Requires: FSL
Requires: Matlab, with the following toolboxes:
signal
Make sure that $FSLDIR/etc/matlab is in your MATLAB path.
Setup: edit the "fix" script so that FIXDIR points at this directory
and so that MATLAB is set to your command line call to matlab
Setup for using inside HCP pipelines (only needed for HCP):
- Edit the "hcp_fix" script in the same way as "fix" described above
- Edit fix_3_clean.m so that variables CIFTI and WBC point to connectome-related files
-------------------------------------------------------
USAGE
See usage instructions at the FSL Wiki:
http://www.fmrib.ox.ac.uk/fsl/FIX
-------------------------------------------------------
function features = featurearfull(timeSeries, TR, nOrder)
% this function calculates the AR-based features for the melodic time series
if(nargin<3), nOrder = 1; end
%timeSeries = iddata(timeSeries,[], TR);
%[m, ~] = ar(timeSeries, nOrder);
%features = [m.NoiseVariance; m.ParameterVector];
[grotAR,grotV] = aryule(timeSeries, nOrder);
features = [grotV grotAR(2:end)];
end
function features = featurearvswgn(timeSeries, TR)
% this function calculates the AR-based features for the melodic time series.
%timeSeries = iddata(timeSeries,[], TR);
%for i = 1:6
% [m, ~] = ar(timeSeries, i);
% mvec(i) = m.NoiseVariance;
%end
for i = 1:6
[grotAR,grotV] = aryule(timeSeries, i);
mvec(i) = grotV;
end
features = polyfit(1:6, mvec, 1);
end
function feat = featureclusterdist(strFold, xSiz)
% this feature tries to quantify the number of clusters in a component
global icThreshold;
call_fsl(sprintf('cluster --in=%sfix/dummyabs --thresh=%f > %sfix/dummy1', strFold, icThreshold, strFold));
counter = 0;
fid = fopen(sprintf('%sfix/dummy1', strFold));
tline = fgetl(fid);
while 1, counter = counter + 1;
tline = fgetl(fid);
if ~ischar(tline), break, end;
C = textscan(tline, '%d %d %f %f %f %f %f %f %f');
drow(counter,:) = [counter C{1} C{2} ];
end;
fclose(fid);
delete(sprintf('%sfix/dummy1*', strFold));
% get all the cluster size
if counter>1
X = double(drow(:,3));
X = X(X>4);
if(length(X)<3), X = [X; zeros(3-length(X),1)]; end
Y = sort(X, 'descend')*xSiz(1)*xSiz(2)*xSiz(3);
feat = [length(X) mean(X)-median(X) max(X) var(X) skewness(X) kurtosis(X) Y(1) Y(2) Y(3)];
else
feat = [0 0 0 0 0 1.5 0 0 0];
end
end
function feat = featureedgemasks(strFold)
% compare high zstat voxels against original EPI mean image intensities
% (in some artefacts the mean FMRI image is dark)
% Similar to above - compare zstat against original EPI *edge* image
% and/or the original PCA residual image
global icThreshold;
feat = [];
for i = 1:5
imgMaskEdge = read_avw(sprintf('%sedge%d',[strFold 'fix/'],i));
imgIc = abs(read_avw(sprintf('%sfix/dummy', strFold)))>icThreshold;
whatPercentOnEdgeB = sum(imgIc(:).*imgMaskEdge(:))/(eps+sum(imgIc(:)));
whatPercentOfEdge = sum(imgIc(:).*imgMaskEdge(:))/(eps+sum(imgMaskEdge(:)));
imgIc = read_avw(sprintf('%sfix/dummyabs', strFold));
wch = imgIc(:)>icThreshold;
whatPercentOnEdgeC = sum(imgIc(wch).*imgMaskEdge(wch))/(eps+sum(imgIc(wch)));
feat = [feat whatPercentOnEdgeB whatPercentOnEdgeC whatPercentOfEdge]; %#ok<AGROW>
end
return
function features = featureentropy(S)
[X, ~] = hist(S, linspace(-2,2,floor(sqrt(length(S)))));
X = X/sum(X);
features(1) = -sum( X.*log(X+eps) );
%features(2) = mean(S.^3)^2/12-kurtosis(S)^2/48; % old broken estimate of negentropy
%features(2) = mean(S.^3)^2/12+(kurtosis(S)-3)^2/48; % corrected but crap estimate of negentropy
S=S(:)/(std(S(:))+eps); features(2)=sum( (mean(-exp(-S.^2/2))+0.71).^2 ); % better-conditioned estimate of negentropy
end
function feat = featureeverynthvarianve(strFold)
% all-of-a-single-slice and none-of-others or every odd (or even) slices -
% effects; the above also tend to show up as VERY sparse or oddly
% nonGaussian in time domain (this is like F4 in Podrack 08)
global icThreshold;
% %
img = read_avw(sprintf('%sfix/dummyabs', strFold));
for i = 1:size(img,3),
imgS = squeeze(img(:,:,i));
varSlc(i) = sum((imgS(:).^2));
end
feat(1) = abs(sum(varSlc(1:2:size(img,3)))-sum(varSlc(2:2:size(img,3))))/(eps+sum(varSlc));
%
AA = [1:4:size(img,3) 2:4:size(img,3)];
BB = [3:4:size(img,3) 4:4:size(img,3)];
feat(2) = abs(sum(varSlc(AA))-sum(varSlc(BB)))/(eps+sum(varSlc));
% %
img = img.*(img>icThreshold);
for i = 1:size(img,3),
imgS = squeeze(img(:,:,i));
varSlc(i) = var(imgS(:));
end
feat(3) = abs(sum(varSlc(1:2:size(img,3)))-sum(varSlc(2:2:size(img,3))))/(eps+sum(varSlc));
%
AA = [1:4:size(img,3) 2:4:size(img,3)];
BB = [3:4:size(img,3) 4:4:size(img,3)];
feat(4) = abs(sum(varSlc(AA))-sum(varSlc(BB)))/(eps+sum(varSlc));
end
function features = featurefftcoarse(timeSeries, TR, freqThreshold)
% In this function, given the time series and its sampling rate (e.g,
% TR), we calculate the ratio of its high-frequency energy to low-frequency
% energy
L = length(timeSeries); % Length of signal
nFFT = 2^nextpow2(L); % Next power of 2 from length of S
Y = abs(fft(timeSeries, nFFT)/L);
f = ((1/TR)/2)*linspace(0,1,nFFT/2);
Y = Y(1:length(f));
features = sum(Y(f>=freqThreshold))/(eps+sum(Y(f<freqThreshold)));
% features = sum(Y(f>=freqThreshold))/(sum(Y)+eps);
end
function features = featurefftfiner(timeSeries, TR, threshArray)
% In this function, given the time series and its sampling rate (e.g,
% TR), we calculate the ratio of its high-frequency energy to low-frequency
% energy
L = length(timeSeries); % Length of signal
nFFT = 2^nextpow2(L); % Next power of 2 from length of S
Y = abs(fft(timeSeries, nFFT)/L);
f = ((1/TR)/2)*linspace(0,1,nFFT/2);
Y = Y(1:length(f));
features = zeros(1,length(threshArray));
for i = 1:(length(threshArray)-1)
features(i) = sum(Y(f>threshArray(i) & f<=threshArray(i+1)));
end
features(i+1) = sum(Y(f>threshArray(i+1)));
features = features/(eps+sum(features));
end
function features = featurefftfinerwrtnull(timeSeries, TR, threshArray)
% % Concolve single-gamma HRF w/WGN and compare its spectrum w/actual RSN
% do it for real
L = length(timeSeries); % Length of signal
nFFT = 2^nextpow2(L); % Next power of 2 from length of S
Y = abs(fft(timeSeries, nFFT)/L);
f = ((1/TR)/2)*linspace(0,1,nFFT/2);
Y = Y(1:length(f));
features = zeros(1,length(threshArray));
for i = 1:(length(threshArray)-1)
features(i) = sum(Y(f>threshArray(i) & f<=threshArray(i+1)));
end
features(i+1) = sum(Y(f>threshArray(i+1)));
featuresReal = features;
% do it under the H0
delay = 6/TR;
sigma = delay/2;
hrf = pdf('gamma', 0:(size(timeSeries,1)-1), (delay/sigma)^2, delay/(sigma^2));
featuresNull = features*0;
for i = 1:100
y = randn(length(timeSeries),1);
y = conv(y,hrf);
y = y(1:length(timeSeries));
y = y/std(y);
yFft = abs(fft(y, nFFT)/L);
yFft = yFft(1:length(f));
for j = 1:(length(threshArray)-1)
features(j) = sum(yFft(f>threshArray(j) & f<=threshArray(j+1)));
end
features(j+1) = sum(yFft(f>threshArray(j+1)));
featuresNull = featuresNull + features;
end
featuresNull = featuresNull/100;
features = (featuresReal-featuresNull).^2./(eps+(featuresReal.^2));
features = [features sum(features)];
end
function features = featuremasktscorrandoverlap(strFold, timeSeries, imgSeg, timeSeriesCsf, timeSeriesGm, timeSeriesWm)
% In this function, given the time series, the mean CSF time-series is
% estimated and then its correlation with the given ts is calc'd
% Next, the overlap of spatial map w/each and everyone of these masks is
% estimated
global icThreshold;
features = regress(timeSeries/(eps+std(timeSeries)), [ones(size(timeSeries)) timeSeriesWm/(eps+std(timeSeriesWm)), timeSeriesGm/(eps+std(timeSeriesGm)), timeSeriesCsf/(eps+std(timeSeriesCsf))]);
features = abs(features(2:end));
imgIc = read_avw(sprintf('%sfix/dummyabs', strFold));
denom = eps+sum(imgIc(:));
features(4) = sum(imgIc(:).*(imgSeg(:)==1))/denom;
features(5) = sum(imgIc(:).*(imgSeg(:)==2))/denom;
features(6) = sum(imgIc(:).*(imgSeg(:)==3))/denom;
imgIc = imgIc>icThreshold;
denom = eps + sum(imgIc(:));
features(7) = sum(imgIc(:).*(imgSeg(:)==1))/denom;
features(8) = sum(imgIc(:).*(imgSeg(:)==2))/denom;
features(9) = sum(imgIc(:).*(imgSeg(:)==3))/denom;
imgIc = read_avw(sprintf('%sfix/dummyabs', strFold));
wch = imgIc>icThreshold;
denom = eps+sum(imgIc(wch));
features(10) = sum(imgIc(wch).*(imgSeg(wch)==1))/denom;
features(11) = sum(imgIc(wch).*(imgSeg(wch)==2))/denom;
features(12) = sum(imgIc(wch).*(imgSeg(wch)==3))/denom;
end
function feat = featuremaxtfce(strFold)
% test for "clusterness" / "compactness" - e.g. feed through TFCE and
% take the max (maybe pre-TFCE normalise e.g. by max(abs(z)))
global icThreshold;
% old version
% system(sprintf('fslmaths %sfix/dummy -abs -div `fslstats %sfix/dummy -a -P 99` -tfce 2 .5 6 %sfix/dummy2', strFold, strFold, strFold));
% new version
call_fsl(sprintf('fslmaths %sfix/dummy -abs -div `fslstats %sfix/dummy -S` -tfce 2 .5 6 %sfix/dummy2', strFold, strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummy2 -P 100 > %sfix/dummy2.txt', strFold, strFold));
feat(1) = load(sprintf('%sfix/dummy2.txt', strFold));
delete(sprintf('%sfix/dummy2*', strFold));
call_fsl(sprintf('fslmaths %sfix/dummy -abs -tfce 2 .5 6 %sfix/dummy2', strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummy2 -P 100 > %sfix/dummy2.txt', strFold, strFold));
feat(2) = load(sprintf('%sfix/dummy2.txt', strFold));
delete(sprintf('%sfix/dummy2*', strFold));
% system(sprintf('fslmaths %sfix/dummy -abs -thr 3 -tfce 2 .5 6 %sfix/dummy2', strFold, strFold));
call_fsl(sprintf('fslmaths %sfix/dummy -abs -thr %f -tfce 2 .5 6 %sfix/dummy2', strFold, icThreshold, strFold));
call_fsl(sprintf('fslstats %sfix/dummy2 -P 100 > %sfix/dummy2.txt', strFold, strFold));
feat(3) = load(sprintf('%sfix/dummy2.txt', strFold));
delete(sprintf('%sfix/dummy2*', strFold));
end
function features = featuremotioncorrelation(timeSeries, strFold)
% calculate the time series' correlation w/motion
thecwd=pwd;
cd(strFold);
[grot,TR]=call_fsl('fslval filtered_func_data pixdim4'); TR=str2num(TR);
[grot,hp]=system('[ -f design.fsf ] && [ _`grep fmri\(temphp_yn\) design.fsf | awk ''{print $3}''` = _1 ] && grep fmri\(paradigm_hp\) design.fsf | awk ''{print $3}''');
hp=str2num(hp);
if grot==1
hp=0;
end
confounds = functionmotionconfounds(TR,hp);
cd(thecwd);
TSS=confounds;
%C1 = corrcoef([timeSeries TSS ]);
%C2 = corrcoef([timeSeries gradient(TSS')']);
%features(1) = max(abs(C1(1,2:end)));
%features(2) = max(abs(C2(1,2:end)));
%features(3) = max(features);
C1 = corrcoef([timeSeries TSS ]); % TSS already contains the derivatives
C1 = abs(C1(1,2:end));
features(1) = max(C1(1:6));
features(2) = max(C1(7:end));
features(3) = max(features);
regModel = regress(timeSeries/(eps+std(timeSeries)), [ones(size(timeSeries)) TSS]);
regModel = sort(abs(regModel(2:end)));
features(4:6) = [regModel(end-1:end); mean(regModel(floor(end/2):end))];
end
function feat = featurenegativevspositive(strFold)
% negatives: if you have a lot of negatives (similar number and
% strength to positives) that's probably an artefact AND general stats
% mean entropy (of nonzero voxels)
global icThreshold;
% system(sprintf('a=`fslstats %sfix/dummy -E`; echo $a>%sfix/dummy1.txt', strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummy -E > %sfix/dummy1.txt', strFold, strFold));
feat(1) = load(sprintf('%sfix/dummy1.txt', strFold));
delete(sprintf('%sfix/dummy1*', strFold));
% system(sprintf('a=`fslstats %sfix/dummyabs -E`; echo $a>%sfix/dummy1.txt', strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummyabs -E > %sfix/dummy1.txt', strFold, strFold));
feat(2) = load(sprintf('%sfix/dummy1.txt', strFold));
delete(sprintf('%sfix/dummy1*', strFold));
% mean (for nonzero voxels)
% system(sprintf('a=`fslstats %sfix/dummy -M`; echo $a>%sfix/dummy1.txt', strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummy -M > %sfix/dummy1.txt', strFold, strFold));
a = load(sprintf('%sfix/dummy1.txt', strFold));
delete(sprintf('%sfix/dummy1*', strFold));
% system(sprintf('a=`fslstats %sfix/dummy -S`; echo $a>%sfix/dummy1.txt', strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummy -S > %sfix/dummy1.txt', strFold, strFold));
aS = load(sprintf('%sfix/dummy1.txt', strFold));
delete(sprintf('%sfix/dummy1*', strFold));
% mean (for abs of nonzero voxels)
% system(sprintf('a=`fslstats %sfix/dummyabs -M`; echo $a>%sfix/dummy1.txt', strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummyabs -M > %sfix/dummy1.txt', strFold, strFold));
b = load(sprintf('%sfix/dummy1.txt', strFold));
delete(sprintf('%sfix/dummy1*', strFold));
% system(sprintf('a=`fslstats %sfix/dummyabs -S`; echo $a>%sfix/dummy1.txt', strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummyabs -S > %sfix/dummy1.txt', strFold, strFold));
bS = load(sprintf('%sfix/dummy1.txt', strFold));
delete(sprintf('%sfix/dummy1*', strFold));
feat(3) = abs(a)/(eps+abs(b));
feat(4) = abs(a*bS)/abs(eps+aS*b);
% diff (for abs of pos and neg voxels)
img = read_avw(sprintf('%sfix/dummy', strFold));
feat(5) = abs(sum(img(:)>0)-sum(img(:)<0))/(eps+sum(img(:)>0));
feat(6) = abs(sum(img(:)>icThreshold)-sum(img(:)<-icThreshold))/(eps+sum(img(:)>icThreshold));
end
function features = featureoupjk(S, deltat)
% The Ornstein Uhlenbeck process is widely used for modelling a
% mean-reverting process. The process S is modelled as
% dS = \lambda (\mu-S)dT + \sigma dW
% Where
% W is a Brownian- Motion, so dW~N(0, dt),
% \lambda meaures the speed of mean reversion
% \mu is the 'long run mean', to which the process tends to revert.
% \sigma, as usual, is a measure of the process volatility
% Since the basic ML has a bias (resulting in frequent estimates of
% lambda which are much too high), we perform a 'jackknife' operation
% to reduce the bias.
% Regressions prefer row vectors to column vectors, so rearrange if
% necessary.
if (size(S,2) > size(S,1)), S = S'; end
m = 2; % Number of partitions
partlength = floor(length(S)/m);
Spart = zeros(m,partlength);
for i=1:1:m
Spart(i,:) = S(partlength*(i-1)+1:partlength*i);
end
[muT, sigmaT, lambdaT ] = functionoupmle(S, deltat);
% Calculate the individual partitions.
mupart = zeros(m,1);
sigmapart = zeros(m,1);
lambdapart= zeros(m,1);
for i=1:1:m
[mupart(i), sigmapart(i), lambdapart(i)] = functionoupmle(Spart(i,:), deltat);
end
% Now the jacknife calculation.
lambda = (m/(m-1))*lambdaT - (sum(lambdapart))/(m^2-m);
% mu and sigma are not biased, so there's no real need for the jackknife.
% But we do it anyway for demonstration purposes.
mu = (m/(m-1))*muT - (sum(mupart ))/(m^2-m);
sigma = (m/(m-1))*sigmaT - (sum(sigmapart ))/(m^2-m);
features = [sigma lambda];
end
function feat = featuresagmasks(strFold)
% look at intersections of ICA spatial map with major vein mask images
global icThreshold;
imgIc2 = read_avw(sprintf('%sfix/dummyabs', strFold));
imgIc = imgIc2>icThreshold;
wch = imgIc2(:)>icThreshold;
feat = [];
for i = 0:3
imgMaskEdge = read_avw(sprintf('%sfix/std1mm2exfunc%d',strFold,i));
whatPercentOnEdgeB = sum(imgIc(:).*imgMaskEdge(:))/max(1,sum(imgIc(:)));
whatPercentOfEdge = sum(imgIc(:).*imgMaskEdge(:))/max(1,sum(imgMaskEdge(:)));
whatPercentOnEdgeC = sum(imgIc2(wch).*imgMaskEdge(wch))/max(1,sum(imgIc2(wch)));
grot1=corrcoef([imgMaskEdge(:) imgIc(:) imgIc2(:)]); grot1(isnan(grot1))=0;
grot2=corrcoef([imgMaskEdge(wch) imgIc2(wch)]); grot2(isnan(grot2))=0;
feat = [feat whatPercentOnEdgeB whatPercentOnEdgeC whatPercentOfEdge grot1(1,2:3) grot2(1,2)];
imgMaskEdge = read_avw(sprintf('%sfix/std1mm2exfunc%ddil',strFold,i));
whatPercentOnEdgeB = sum(imgIc(:).*imgMaskEdge(:))/max(1,sum(imgIc(:)));
whatPercentOfEdge = sum(imgIc(:).*imgMaskEdge(:))/max(1,sum(imgMaskEdge(:)));
whatPercentOnEdgeC = sum(imgIc2(wch).*imgMaskEdge(wch))/max(1,sum(imgIc2(wch)));
grot1=corrcoef([imgMaskEdge(:) imgIc(:) imgIc2(:)]); grot1(isnan(grot1))=0;
grot2=corrcoef([imgMaskEdge(wch) imgIc2(wch)]); grot2(isnan(grot2))=0;
feat = [feat whatPercentOnEdgeB whatPercentOnEdgeC whatPercentOfEdge grot1(1,2:3) grot2(1,2)];
imgMaskEdge = read_avw(sprintf('%sfix/std1mm2exfunc%ddil2',strFold,i));
whatPercentOnEdgeB = sum(imgIc(:).*imgMaskEdge(:))/max(1,sum(imgIc(:)));
whatPercentOfEdge = sum(imgIc(:).*imgMaskEdge(:))/max(1,sum(imgMaskEdge(:)));
whatPercentOnEdgeC = sum(imgIc2(wch).*imgMaskEdge(wch))/max(1,sum(imgIc2(wch)));
grot1=corrcoef([imgMaskEdge(:) imgIc(:) imgIc2(:)]); grot1(isnan(grot1))=0;
grot2=corrcoef([imgMaskEdge(wch) imgIc2(wch)]); grot2(isnan(grot2))=0;
feat = [feat whatPercentOnEdgeB whatPercentOnEdgeC whatPercentOfEdge grot1(1,2:3) grot2(1,2)];
end
return
function feat = featureslicewisestats(strFold)
% all-of-a-single-slice and none-of-others or every odd (or even) slices -
% effects; the above also tend to show up as VERY sparse or oddly
% nonGaussian in time domain
global icThreshold;
img = read_avw(sprintf('%sfix/dummyabs', strFold));
sumTot = sum(img(:).^2)+eps;
for i = 1:size(img,3)
imgS = squeeze(img(:,:,i));
sumSlc(i) = sum(imgS(:).^2)/sumTot;
end
feat(1) = sum(sumSlc>.15);
feat(2) = max(sumSlc);
sumTot = sum(img(:)>icThreshold)+eps;
for i = 1:size(img,3)
imgS = squeeze(img(:,:,i));
sumSlc(i) = sum(imgS(:)>icThreshold)/sumTot;
end
feat(3) = sum(sumSlc>.15);
feat(4) = max(sumSlc);
end
function feat = featuresmoothest(strFold)
% closely related: spatial smoothness (resel size etc) of zstat map
% spatial entropy (histogram) / parameters from histogram mixture-model fit
% smoothness
msk = sprintf('%smask',strFold);
call_fsl(sprintf('echo `smoothest -z %sfix/dummy -m %s` > %sfix/dummy1.txt', strFold, msk, strFold));
fid = fopen(sprintf('%sfix/dummy1.txt', strFold));
tline = fgetl(fid);
C = textscan(tline, '%s %f %s %f %s %f');
fclose(fid);
% DLH, VOL, RESEL
feat = C{6};
delete(sprintf('%sfix/dummy1*', strFold));
end
function feat = featurestripe(strFold)
call_fsl(sprintf('fslmaths %sfix/dummy -s 2 -abs %sfix/dummy01', strFold, strFold));
call_fsl(sprintf('fslmaths %sfix/dummy -abs -s 2 -sub %sfix/dummy01 -abs %sfix/dummy01', strFold, strFold, strFold));
call_fsl(sprintf('fslstats %sfix/dummy01 -P 95 > %sfix/dummy01', strFold, strFold));
feat = dlmread(sprintf('%sfix/dummy01', strFold));
delete(sprintf('%sfix/dummy0*', strFold));
end
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment