Creating Paired Differences

If you analyze 2-time point longitudinal data, you will eventually observe that it is easier to create and analyze difference data (e.g. data from time 2 minus data from time 1). However, if you’re not a scripting guru this might be an annoying task.

This PairDiff.m script will take an even number of images and produce a set of paired difference images, even minus odd.

function PairDiff(Imgs,BaseNm)
% PairDiff(Imgs,BaseNm)
% Imgs   - Matrix of (even-number of) filenames
% BaseNm - Basename for difference images
%
% Create pairwise difference for a set of images (img2-img1, img4-img3, etc),
% named according to BaseNm (BaseNm_0001, BaseNm_0002, BaseNm_0003, etc).
%______________________________________________________________________
% $Id: PairDiff.m,v 1.2 2012/02/16 21:36:41 nichols Exp $

if nargin<1 | isempty(Imgs)
  Imgs = spm_select(Inf,'image','Select n*2 images');
end
if nargin<2 | isempty(BaseNm)
  BaseNm = spm_input('Enter difference basename','+1','s','Diff');
end
V = spm_vol(Imgs);
n = length(V);
if rem(n,2)~=0
  error('Must specify an even number of images')
end

V1 = V(1);
V1  = rmfield(V1,{'fname','descrip','n','private'});
V1.dt = [spm_type('float32') spm_platform('bigend')];

for i=1:n/2
  V1.fname = sprintf('%s_%04d.img',BaseNm,i);
  V1.descrip = sprintf('Difference %d',i);
  Vo(i) = spm_create_vol(V1);
end

%
% Do the differencing
%

fprintf('Creating differenced data ')

for i=1:2:n

  img1 = spm_read_vols(V(i));
  img2 = spm_read_vols(V(i+1));

  img = img2-img1;

  spm_write_vol(Vo((i+1)/2),img);

end

After I created this I realized that Ged Ridgway also has a similar script, make_diffs.m , that takes two lists of images (baseline, followup) and does the same thing, though with perhaps more intuitive filenames.

Flame without 1st level directories – copes only.

Follow-up to Flame without 1st level directories from the NISOx blog (formerly Neuroimaging Statistics Tips & Tools)

My earlier post described how to get the basic FEAT second-level results when you just have a stack of COPE and VARCOPE images. In fact, you don’t even need the VARCOPE images.

If you just have the COPE images, in the flameo call omit --vc=4dvarcope (which you don’t have) and replace the --runmode=flame1 with --runmode=ols and you should be all set.

Save As HTML in Powerpoint 2011 for Mac

Follow-up to Extracting Images from Powerpoint Presentations for (e.g.) LaTeX presentations from the NISOx blog (formerly Neuroimaging Statistics Tips & Tools)

[NOTE: Powerpoint can now do this directly via “Export as…”, and then selecting HTML. Leaving this entry just for reference. -TN 2017/04/27]

Arrrgh. In Powerpoint 2011 for Mac there is no “HTML” option in the “Save as…” dialog.

Fortunately, there is a Microsoft Support note that gives a work around. For convenience (& my memory) here it is distilled:

  1. Press ALT+F11. (Depending on your keyboard settings, you might have to hold down “fn” to get F11.) This will open a “Project” window.
  2. Press CMD+CTRL+G to open the “Immediate” window. (The support page leaves out CMD, the ‘Apple’ key.)
  3. In the “Immediate” window paste in the following:
    ActivePresentation.SaveAs "/tmp/name.htm", ppSaveAsHTML, msoFalse
    Press ENTER.

This will give you a name.htm file and a name_files directory. Of course, you can replace ”/tmp/myfile.htm” with the desired path and filename. However, I find it saves keystrokes to save it into /tmp, and then move the file to wherever you really want it.

While you’ll see the images in the *_files/ folder, disappointingly, I don’t find the same original source images like I have seen before. Again, Arrrgh.

-Tom

SPM8 Gem 1: Zero NaN’s with the zeronan.m script

Follow-up to SPM99 Gem 3: NaNing zero values from the NISOx blog (formerly Neuroimaging Statistics Tips & Tools)

This was the topic of SPM99 Gem 3, converting NaN’s to zeros. For SPM8, see the following script zeronan.m that will zero NaN’s for you.

-Tom

function ofNm = zeronan(ifNm,val)
% FORMAT ofNm = zeronan(ifNm,val)
% ifNm  - Input filename(s)
% val   - Value to set NaN's to (defaults to zero)
%
% Output:
% ofNm  - Cell array of output filenames.
%
%
% Images have NaN's replaced with zeros, and new images, prefixed with a
% 'z', are created.
%
%________________________________________________________________________
% Based on zeronan.m,v 1.3 2005/10/26 21:58:55 nichols Exp 
% Thomas Nichols, 1 April 2011

if nargin<2, val = 0; end
if nargin0'); end

if ~iscell(ifNm)
  ifNm = cellstr(ifNm)';
else
  ifNm = ifNm(:)';
end

OtfNm = {};

for fNm = ifNm

  fNm = fNm{:};

  OfNm = ['z' fNm];
  [pth,nm,xt,vol] = spm_fileparts(fNm);
  OfNm = fullfile(pth,['z' nm xt]);

  % Code snippet from John Ashburner...
  VI       = spm_vol(fNm);
  VO       = VI;
  VO.fname = OfNm;
  VO       = spm_create_vol(VO);
  for i=1:VI.dim(3),
    img      = spm_slice_vol(VI,spm_matrix([0 0 i]),VI.dim(1:2),0);
    tmp      = find(isnan(img));
    img(tmp) = val;
    VO       = spm_write_plane(VO,img,i);
  end;

  OtfNm = {OtfNm{:} OfNm};

end

if nargout>0
  ofNm = OtfNm;
end

Extracting Images from Powerpoint Presentations for (e.g.) LaTeX presentations

Have you ever needed to create a LaTeX presentation but needed figures that only lived in a PowerPoint file? Or have you found that the only copy of a crucial figure you need lives inside a PowerPoint presentation, and when you copy it out of PowerPoint into another application it looks horrific?

This is a cool trick to extract the original (or close to original) image file from the presentation.

Open the presentation, delete all the slides you don’t need. The option is different with different versions of Powerpoint; look for “File -> Export as…” or ” “File -> Save as Web Page…”. You will get a .html (or .htm) file, and a _files directory. Inside the _files directory you will find a bunch of image files. It might take some poking around if there are lots of images on the slide, but you should find a close to original version of the image pasted into the slide.

Orthogonal polynomials with intuitive scale

When creating regressors (or contrasts) for a parametric experimental design, it is often useful to model (or test) linear, quadratic and cubic effects separately. However, if you just use the usual x, x*x, x*x*x approach to create predictors, you wil create highly correlated variables. Specifically, the quadratic term will include some of the linear effect, and vice-versa. The solution is to use orthogonal polynomials.

The very simplest way to do this is with Matlab’s orth function. E.g. for 10-level parametric regressor

n=10;
x=[1:n]';
One=ones(n,1);
X0=[One x x.^2 x.^3];
X=orth(X0);
plot(X)

But annoyingly, this standardizes the columns and produces a fairly uninterprtable basis. In particular, while a linear trend can be represented by a linear combination of those columns, it is not in any one single column.

Instead, we can carefully build up a basis that is mean zero, and uses successive othogonalization to ensure the interpretability of each column.

n=10;
x=[1:n]';
One=ones(n,1);
X0=[x x.^2 x.^3];
X=One;
for i=1:size(X0,2)
  X=[X X0(:,i)-X*pinv(X)*X0(:,i)];
end
plot(X)

Now we have three columns readily identifiable as intercept, linear, quadratic and cubic terms; after the first column, they are mean zero (so that they do not carry any of the grand mean), and all are orthogonal, as can be checked with corrcoef(X).

There are of course specialized functions to create orthogonal polynomials, but for this purpose, this code snippet is all you need.

Flame without 1st level directories

FSL’s fMRI modelling tool, FEAT, is a great tool and is easy to use if you use FSL from start to finish, in the recommended fashion. However, sometimes you might want to get FEAT group results when you don’t have all of the individual 1st level FEAT directories. This tip is for the case when you have just the 4D COPE image and the 4D VARCOPE images (all subjects, of course, already in MNI space), and you just want to do a one-sample t-test.

In the code below, I’m assuming that the following files are around:

  • 4dcope 4D COPE image
  • 4dvarcope 4D VARCOPE image
  • design.mat, design.con and design.con for a 1-sample t-test; use Glm_gui if you need help creating these.
  • mask mask image
  • example_func an example image for ‘underlay’ in the figures
$FSLDIR/bin/flameo --cope=4dcope --vc=4dvarcope --mask=mask --ld=stats --dm=design.mat --cs=design.grp --tc=design.con --runmode=flame1
echo $($FSLDIR/bin/fslnvols 4dcope) - 1 | bc -l  > stats/dof
/bin/rm -f stats/zem* stats/zols* stats/mask*
$FSLDIR/bin/smoothest -d $(cat stats/dof) -m mask -r stats/res4d > stats/smoothness
rm -f stats/res4d*
awk '/VOLUME/ {print $2}' stats/smoothness > thresh_zstat1.vol
awk '/DLH/ {print $2}' stats/smoothness > thresh_zstat1.dlh
$FSLDIR/bin/fslmaths stats/zstat1 -mas mask thresh_zstat1
$FSLDIR/bin/cluster -i thresh_zstat1 -c stats/cope1 -t 2.3 -p 0.05 -d $(cat thresh_zstat1.dlh) --volume=$(cat thresh_zstat1.vol) --othresh=thresh_zstat1 -o cluster_mask_zstat1 --connectivity=26 --mm --olmax=lmax_zstat1_tal.txt > cluster_zstat1_std.txt
$FSLDIR/bin/cluster2html . cluster_zstat1 -std
MinMax=$($FSLDIR/bin/fslstats thresh_zstat1 -l 0.0001 -R)
$FSLDIR/bin/overlay 1 0 example_func -a thresh_zstat1 $MinMax rendered_thresh_zstat1
$FSLDIR/bin/slicer rendered_thresh_zstat1 -S 2 750 rendered_thresh_zstat1.png
/bin/cp $FSLDIR/etc/luts/ramp.gif .ramp.gif

You don’t get the full, pretty output from FEAT, but you do get a pretty PNG overlay, and the HTML table of significant clusters. Generalizing to more complicated group-level models is simply a matter of changing the design.mat and design.con files, and adjusting the degrees-of-freedom calculation on the second line.

This is really just a quick hack. I hope it is of some use.

Quick addition: The code above gives you corrected cluster-wise inference. To get corrected voxel-wise inference, replace the line with the ”$FSLDIR/bin/cluster” command with these four lines:

RESELcount=$(awk '$1~/VOLUME/{v=$2};$1~/RESELS/{r=$2};END{printf("%g",1.0*v/r)}' stats/smoothness)
FWEthresh=$(ptoz 0.05 -g $RESELcount)
$FSLDIR/bin/fslmaths thresh_zstat1 -thr $FWEthresh thresh_zstat1
$FSLDIR/bin/cluster -i thresh_zstat1 -c stats/cope1 -t $FWEthresh --othresh=thresh_zstat1 -o cluster_mask_zstat1 --connectivity=26 --mm --olmax=lmax_zstat1_tal.txt > cluster_zstat1_std.txt

Or, if you want to use an uncorrected threshold, use this snippet instead (replacing 2.3 with your favorite threshold).

UnCorrThresh=2.3
$FSLDIR/bin/cluster -i thresh_zstat1 -c stats/cope1 -t $UnCorrThresh --othresh=thresh_zstat1 -o cluster_mask_zstat1 --connectivity=26 --mm --olmax=lmax_zstat1_tal.txt > cluster_zstat1_std.txt