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).

$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

Make Word act like Emacs

Follow-up to Mac Goodness from the NISOx blog (formerly Neuroimaging Statistics Tips & Tools)

I’ve gone on about how great it is to change keyboard shortcuts with System Preferences -> Keyboards, but not all Word commands are accessible through that facility. In particular, it is possible to make Word respond to Emacs editing commands!

In Word, go to Tools -> Customize Keyboard… . Under Categories, select “All Commands” to show everything. In the Commands box find a command, and then add the corresponding key as follows:

  • Command StartOfLine – Shortcut Control+A
  • Command EndOfLine – Shortcut Control+E
  • Command EditCut – Shortcut Control+W
  • Command EditPaste – Shortcut Control+Y
  • Command EditFind – Shortcut Control+S

Mac Tweaks

This is just a rambling list of the tweaks I’ve made to my Mac, mostly to help me remember them if my Mac ever dies, but others might find them useful too.

Disable Spotlight arithmetic

Due to editorial work, I’m constantly using Spotlight to find filenames consisting of numbers like 2010-0021, for which Spotlight helpfully returns the top answer as 1989 (I prefer to use bc on the command line for such arithmetic). Since this is the top answer, if you press return it opens the Calculator program instead of the file you wanted. (ug).

Instead, you can just tell Spotlight not to do any such arithmetic by using this command in the terminal:

defaults write com.apple.spotlight CalculationEnabled NO

You need to kill Spotlight and have it restart to take effect. (While you’re supposed to be able to sudo killall ... to get it to restart, it didn’t work for me, and I just logged out and logged back in).

Force Preview to print at 100% (no scaling) by default

By default, when you print a document in Preview it will ‘scale to fit’. While this seems like a sensible thing, it bases the default scaling on the full size of what you’re printing and the printable area of the paper used. While this is sensible for a JPEG photo, it is idiotic for PDF document.

Finally I discovered that this behavior is easy to change. Just issue this command in a Terminal window

defaults write com.apple.Preview PVImagePrintingScaleMode 0

There is only one down side: If you then change “Images per page” to 2 or greater, the scaling is 100% and the image will be cropped. You just need to select “Scale to fit” (the previous default) and you’ll get the anticipated (appropriately reduced) result.

Prevent tar from saving extended information (resource fork) in ._ files

Set this in your environment:

Mac Goodness

Essential tips to increase my Mac love

Mac is the unification of Unix and an amazingly elegant user interface. Still, there are some keyboard-centric tricks that are not well documented. These are the Mac tips/tricks that bring me happiness.

  1. Change directory from file selection dialog box: Press ~ to pop-up a text box to enter a filepath.
  2. Cycle through an application’s windows with the same shortcut, CMD-~.
  3. I never, never use that yellow minimize button. Instead, I always hide with CMD-h.
  4. If you’re an emacs user, you know the vital importance of the control key. Thus, I’ve used System Preferences -> Keyboard -> Keyboard -> Modifier Keys to remap CAPS LOCK to be an additional CONTROL key. Magic!
    Now if I could only remap the stupid UK-keyboard’s plus/minus/section(!) key to be ESCAPE. Sigh.
  5. Add those missing keyboard short cuts with System Preferences -> Keyboard -> Keyboard Shortcuts. Some of my favorites
    • MS Office apps: CMD-SHIFT-V for Paste Special…
    • MS Office Powerpoint: CTRL-= for Subscript, CTRL-+ for Superscript, though first you need to add these as menu items with View -> Customise Toolbars and Menus; it doesn’t matter which menu you add them to (though Format is the obvious one).
    • MS Office aps: Auto Correct common LaTeX notation into symbols. One time, type the symbol you want (e.g. β), copy it. Under Tools->Autocorrect, enter $beta$ in the “Replace:” field, and paste the symbol into the “With:” field… et voila! The next type you type ”$beta$” it will automagically be replaced with a ”β”. Add other commonly used symbols as needed. (This also works on PC Office apps).
    • Terminal Program: CMD-LEFT and CMD-RIGHT for Select Next Tab, Select Previous Tab. (And, if you haven’t discovered the magic of Tabs in Terminal program, you’re just a CMD-t away from bliss.)
  6. Having trouble with the keyboard short cuts not working? The text of the menu item must be identical to what you enter in the keyboard shortcut dialog; note that an ellipsis (…) is not three periods (…); get an ellipses on the Mac with OPTION-;
  7. Terminal bliss. I almost never use the Finder, as the Terminal can do everything I need and faster. There are, on occasions, though, things that the Finder is better at. On those occasions I just open the directory, like open ., which gives me a finder window for the current directory in the shell, or open ~/Talks/Mytalk for another directory.
  8. Terminal bliss. I find ejecting a volume cumbersome… you have to mouse around, find the volume, click on it, hit CMD-e. Ug. Better, is just to issue the command hdiutil eject /Volume/<VolumeToEject>. Better yet, define the following alias: alias eject hdiutil eject. Then you can intuitively just type eject /Volume/<VolumeToEject>. Magic!
  9. Terminal bliss. I have a “sleepnow” alias so I don’t have to fumble for the Apple->Sleep menu nor remember to eject my backup volume:
    ejectbackup; osascript -e 'tell application "System Events" to sleep'; echo Going to sleep, please wait...
    The command ejectbackup is yet another alias for:
    if (-d '/Volumes/My Backup') hdiutil eject '/Volumes/My Backup'
    which is c-shell specific, though a bash user would want
    if [ -d '/Volumes/My Backup' ] hdiutil eject '/Volumes/My Backup'

    While this ejectbackup bit is optional, it’s handy since it ejects my Time Machine backup volume so that I can yank the USB cable just before it sleeps; otherwise, without it, you’ll wake up the Mac when you pull the cable and get pesky messages telling you to eject the USB properly. (Of course, replace ‘My Backup’ for the name of your backup volume.)

  10. Remove Apple Extended Attributes. MacOS looks like Unix, but one difference is that Apple has created extensions to record arbitrary key-value pair attributes for each file. For example, the way that your Mac knows you downloaded something from the internet is because it has added attributes with keys like “com.apple.quarantine” to flag the origin of the file. To view such attributes you need the xattr program (from e.g. here), as in
    xattr -l file.txt
    To clear one particular attributes use
    xattr -d <key> file.txt
    where <key> is “com.apple.quarantine” or whatever you want to delete

    To wipe out all attributes you have to do it one by one, ug. Here’s a one-liner to do it

    xattr --list-parseable file.txt | grep = | awk -F= '{print $1}' | xargs -I %% xattr -d %% file.txt

    and, specifically, my csh alias

    alias xattr_del_all 'xattr --list-parseable \!:1 | grep = | awk -F= '\''{print $1}'\'' | xargs -I %% xattr -d %% \!:1'

  11. Sparse Bundles for extra security. I find it handy to keep confidential information on sparse bundles protected with a password. To create one use this command

    hdiutil create -size 100m -volname PatientInfo -encryption -type SPARSEBUNDLE -fs HFS+J ~/Project/PatientInfo

    You’ll be prompted for a password, and then you’ll see that this will create a directory “PatientInfo.sparsebundle” in the Project directory. In order to open this directory you’ll have to enter the password.

    If you ever need to increase the size later, use this command

    hdiutil resize -size 200m ~/Project/PatientInfo.sparsebundle

    To later change the password, use

    hdiutil chpass ~/Project/PatientInfo.sparsebundle

  12. Reunite resource fork on un-archiving, OR, how to avoid __MACOS folders. Mac OS files have a data and resource fork, where the resource fork has information metadata like file creation and modification time. Whenever you unzip a .zip file on the command line that was created on the Mac, you’ll find a __MACOS folder gets created which which contains resource for information in separate ._* files, and they can generally be safely deleted. But if you want to reunite the resource fork with each file, use ditto, like ditto -xk <<ARCHIVE.ZIP>> <<TargetDirectory>>, for example
    ditto -xk Archive.zip .
  13. Resize images. Modern digital cameras create gigantic image files, not great for sharing. While there are variety of ways to resize images on the command line (e.g. with ImageMagick) I was amazed to learn that there’s a tool built right in to Mac OS to do this. The sips command is a “scriptable image processing system”, and to resize images just do

    sips -Z 1024 *.JPG

    but note that this resizes in place, i.e. writes over the existing image. Here, 1024 specifies the maximal dimension the resized images have.

SPM5 Gem 6: Corrected cluster size threshold

This is SPM5 version of SPM2 Gem 13.

This is a script that will tell you the corrected cluster size threshold for given cluster-defining threshold: CorrClusTh.m

The usage is pretty self explanatory:

 Find the corrected cluster size threshold for a given alpha
 function [k,Pc] =CorrClusTh(SPM,u,alpha,guess)
 SPM   - SPM data structure
 u     - Cluster defining threshold
         If less than zero, u is taken to be uncorrected P-value
 alpha - FWE-corrected level (defaults to 0.05)
 guess - Set to NaN to use a Newton-Rhapson search (default)
         Or provide a explicit list (e.g. 1:1000) of cluster sizes to
         search over.
         If guess is a (non-NaN) scalar nothing happens, except the the
         corrected P-value of guess is printed. 

 Finds the corrected cluster size (spatial extent) threshold for a given
 cluster defining threshold u and FWE-corrected level alpha. 

To find the 0.05 (default alpha) corrected cluster size threshold for a 0.01 cluster-defining threshold:

>> load SPM
>> CorrClusTh(SPM,0.01)
  For a cluster-defining threshold of 2.4671 the level 0.05 corrected
  cluster size threshold is 7860 and has size (corrected P-value) 0.0499926

Notice that, due to the discreteness of cluster sizes you cannot get an exact 0.05 threshold.

The function uses an automated search which may sometimes fail. If you specify a 4th argument you can manually specify the cluster sizes to search over:

>> CorrClusTh(SPM,0.01,0.05,6000:7000)

  WARNING: Within the range of cluster sizes searched (6000...7000)
  a corrected P-value <= alpha was not found (smallest P: 0.0819589)

  Try increasing the range or an automatic search.

Ooops… bad range.

Lastly, you can use it as a look up for a specific cluster size threshold. For example, how much over the 0.05 level would a cluster size of 7859 be?

>> CorrClusTh(SPM,0.01,0.05,7859)
  For a cluster-defining threshold of 2.4671 a cluster size threshold of
  7859 has corrected P-value 0.050021

Just a pinch!

SPM5 Gem 1: Introduction to SPM5 scripting

Date: Tue, 13 Dec 2005 17:03:13 +0100
From: John Ashburner <john@FIL.ION.UCL.AC.UK>
Subject: Re: [SPM] Where can we find some development materials for SPM ?

>    As you know,we usually need to modify the code of SPM to fit our
> problem.but we can not find some relevant development  tutorials. Would you
> please tell me how to learn the framework of SPM step by step ?
>   In addition, I want to know where I can find the details of the SPM
> structure.

It may be easiest to learn by example.  If you want to develop a new 
user-interface for SPM5, then you would create a file called spm_config_*.m, 
similar to the other spm_config.m files (if you strip out the documentation 
parts, you will see that these are actually quite small).  Your spm_config* 
file can then be added to the toolbox subdirectory and accessed through the 
TOOLS pulldown.

The help button allows you to navigate through the documentation of each of 
the Matlab functions, which you may find useful.

For reading and writing images, you would use these functions.



There is Matlab help on all these functions.  Alternatively, you could use the 
NIFTI routines directly.  There is no documentation on this, but here are a 
few examples of how you can use them:

% Example of creating a simulated .nii file.
dat            = file_array;
dat.fname = 'junk.nii';
dat.dim     = [64 64 32];
dat.dtype  = 'FLOAT64-BE';
dat.offset  = ceil(348/8)*8;

% alternatively:
% dat = file_array( 'junk.nii',dim,dtype,off,scale,inter)


% Create an empty NIFTI structure
N = nifti;

fieldnames(N) % Dump fieldnames

% Creating all the NIFTI header stuff
N.dat = dat;
N.mat = [2 0 0 -110 ; 0 2 0 -110; 0 0 -2 92; 0 0 0 1];
N.mat_intent = 'xxx'
N.mat_intent = 'Scanner';
N.mat0 = N.mat;
N.mat0_intent = 'Aligned';

N.diminfo.slice = 3;
N.diminfo.phase = 2;
N.diminfo.frequency = 2;
N.diminfo.slice_time.code = 'sequential_increasing';
N.diminfo.slice_time.start = 1;
N.diminfo.slice_time.end = 32;
N.diminfo.slice_time.duration = 3/32;

N.intent.code='xxx' ; % dump possibilities
N.intent.code='FTEST'; % or N.intent.code=4;
N.intent.param = [4 8];

N.timing.toffset = 28800;
N.descrip = 'This is a NIFTI-1 file';
N.cal = [0 1];

create(N); % Writes hdr info


[i,j,k] = ndgrid(1:64,1:64,1:32);
dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 30^2))=1;
dat(find((i-32).^2+(j-32).^2+(k*2-32).^2 < 15^2))=2;

% displaying a slice

% get a handle to 'junk.nii';


Best regards,