Hyo Jong wrote:
> Hello,
> A batch file for Spm5 by R. Henson calls spm5spike.
> I am wondering if anyone has the spm5spike tool and share with us.
> It seems very nice to check the outliers.
>
> Thanks in advance
> Hyo Jong Lee
>
--
-------------------------------------------------------
DR RICHARD HENSON
MRC Cognition & Brain Sciences Unit
15 Chaucer Road
Cambridge, CB2 7EF
England
EMAIL: [log in to unmask]
URL: http://www.mrc-cbu.cam.ac.uk/people/rik.henson/personal
TEL +44 (0)1223 355 294 x522
FAX +44 (0)1223 359 062
MOB +44 (0)794 1377 345
-------------------------------------------------------
function [out,out_scans,dout,dout_scans] = spm5spike(P,thr,pflag,sflag,normflag);
if nargin<5
normflag=1; % divide by std over scans
if nargin<4
sflag=0; % save flag
if nargin<3
pflag=1; % print (graph) flag
end
end
end
% Report 0th and 1st order moments of scans (slice by slice)
% Uses raw data and finited differenced data.
% works only for transaxially oriented scans
% Christian Buechel, Oliver Josephs, Rik Henson
Fgraph = spm_figure('FindWin','Graphics');
if(nargin<1)
[P, flag] = spm_select(Inf,'image');
end
if ischar(P)
try
V=spm_vol(P);
catch
error('input list of files or memory mapped list of files')
end
else
V=P;
end
q = size(P,1);
V = spm_vol(P);
if(nargin<2)
thr = spm_input('threshold in sd',1);
end
dim = V(1).dim; % assumes images all same format/dimensions
Mean_result = zeros(dim(3),q);
Min_result = zeros(dim(3),q);
Max_result = zeros(dim(3),q);
Std_result = zeros(dim(3),q);
%-Cycle over planes
%-----------------------------------------------------------------------
spm_progress_bar('Init',100,'Spike',' ');
for i = 1:dim(3)
B = spm_matrix([0 0 i]);
dat = zeros(prod(V(1).dim(1:2)), q);
msk = logical(zeros(prod(V(1).dim(1:2)),1));
for j = 1:q
M = V(j).mat\V(1).mat*B;
d = spm_slice_vol(V(j),M,V(1).dim(1:2),[1,NaN]);
msk = msk | ~finite(d(:));
dat(:,j) = d(:);
end
dat = dat(find(~msk),:);
if ~isempty(dat),
Mean_result(i,:) = mean(dat);
Std_result(i,:) = std(dat);
end;
clear msk dat
spm_progress_bar('Set',i*100/dim(3));
end % (loop over planes)
spm_progress_bar('Clear');
% Added by Rik to stop normit failing when a slice has zero mean
% (eg due to movement out of field of view)
%tmp = var(Mean_result').*var(Std_result');
%valid_slices = find(tmp>eps);
valid_slices = find(var(Mean_result')>eps & var(Std_result')>eps)
invalid = setdiff(valid_slices,1:dim(3));
if ~isempty(invalid)
sprintf('Slices with zero variance:')
disp(invalid)
end
Mean_result = Mean_result(valid_slices,:);
Std_result = Std_result(valid_slices,:);
Meanmin=min(min(Mean_result));
Meanmax=max(max(Mean_result));
Stdmin=min(min(Std_result));
Stdmax=max(max(Std_result));
%Normalize to mean = 0, std = 1
%---------------------
Nmean = normit(Mean_result',normflag);
Nstd = normit(Std_result',normflag);
Ndmean = normit(diff(Mean_result'),normflag);
Ndstd = normit(diff(Std_result'),normflag);
Ex_mean = find(abs(Nmean) > thr);
Ex_std = find(abs(Nstd) > thr);
Ex_dmean = find(abs(Ndmean) > thr);
Ex_dstd = find(abs(Ndstd) > thr);
Mask_mean = zeros(size(Nmean));
Mask_std = zeros(size(Nstd));
Mask_dmean = zeros(size(Ndmean));
Mask_dstd = zeros(size(Ndstd));
Mask_mean(Ex_mean) = ones(size(Ex_mean));
Mask_std(Ex_std) = ones(size(Ex_std));
Mask_dmean(Ex_dmean) = ones(size(Ex_dmean));
Mask_dstd(Ex_dstd) = ones(size(Ex_dstd));
%display results
%--------------------------------------
if(pflag)
%Fgraph=1; % in SPM2?
%Mean
%----
%figure(Fgraph); spm_clf;
figure;
%colormap hot;
%subplot(6,1,1);
%imagesc(Mean_result);
%caxis([-thr thr]);
%caxis([Meanmin Meanmax]);
%colorbar;
%xlabel('scan');
%ylabel('slice');
%title('Mean');
%Std
%---
%subplot(6,1,2);
%imagesc(Std_result);
%caxis([Stdmin Stdmax]);
%caxis([-thr thr]);
%colorbar;
%xlabel('scan');
%ylabel('slice');
%title('Std');
%NMean
%-----
%subplot(6,1,3);
%imagesc(Nmean');
%caxis([-thr thr]);
%colorbar;
%xlabel('scan');
%ylabel('slice');
%title('Normalised Mean');
title(sprintf('Spikes %s', pwd))
hold on
subplot(4,1,1)
plot(Nmean)
ylabel('slice');
title('Normalised Mean');
%NStd
%----
%subplot(6,1,4);
%imagesc(Nstd');
%caxis([-thr thr]);
%colorbar;
%xlabel('scan');
%ylabel('slice');
%title('Normalised Std');
subplot(4,1,2)
plot(Nstd)
ylabel('slice');
title('Normalised Std');
%dMean
%----
%figure(Fgraph);
%colormap hot;
%subplot(6,1,5);
%imagesc(Ndmean');
%caxis([-thr thr]);
%colorbar;
%xlabel('scan');
%ylabel('slice');
%title('Finite Differenced Mean');
subplot(4,1,3)
plot(Ndmean)
ylabel('slice');
title('Finite Differenced Mean');
%dStd
%---
%subplot(6,1,6);
%imagesc(Ndstd');
%caxis([-thr thr]);
%colorbar;
%xlabel('scan');
%ylabel('slice');
%title('Finite Differenced Std');
subplot(4,1,4)
plot(Ndstd)
xlabel('scan');
ylabel('slice');
title('Finite Differenced Std');
end
%-Unmap volumes
%-----------------------------------------------------------------------
%for i = 1:q; spm_unmap(V(:,i)); end
fprintf('\nSuspicious scans and slices on mean (threshold = %1.1f):',thr);
out = []; % 0th order outliers
out_scans = [];
for j = 1:q
new = 1;
for i = 1:length(valid_slices)
if Mask_mean(j,i) | Mask_std(j,i)
out = [out; j valid_slices(i)];
if new
out_scans = [out_scans; P(j,:)];
fprintf(['\nScan: %1.0f [' P(j,:) ']' '\nSlices (Mean/Std): '], j);
new = ~new;
end
fprintf(['%1.0f (%1.0f,%1.0f) '],valid_slices(i),Nmean(j,i),Nstd(j,i));
end
end
end
fprintf('\n\nSuspicious scans (N+1) and slices on differences (threshold = %1.1f):',thr);
% If diff(j) is large then assume scan n+1 is cause.
% If only n+1
dout = []; % 1st order outliers
dout_scans = [];
for j = 1:(q-1)
new = 1;
for i = 1:length(valid_slices)
if Mask_dmean(j,i) | Mask_dstd(j,i)
dout = [dout; j+1 valid_slices(i)];
if new
dout_scans = [dout_scans; P(j+1,:)];
fprintf(['\nScan: %1.0f [' P(j+1,:) ']' '\nSlices (Mean/Std): '], j+1);
new = ~new;
end
fprintf(['%1.0f (%1.0f,%1.0f) '],valid_slices(i),Ndmean(j,i),Ndstd(j,i));
end
end
end
fprintf('\n');
if sflag
save(fullfile(spm_str_manip(P(1,:),'h'),'/spikes.mat'),'out','out_scans','dout','dout_scans')
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function nv=normit(m,normflag)
for v=1:size(m,2)
nv(:,v)=detrend(m(:,v),0);
if normflag
nv(:,v)=nv(:,v)/std(nv(:,v));
end
end
|