NIfTI-1 is a data format proposed by the NifTI (Neuroimaging Informatics Technology Initiative) DFWG (Data Format Working Group) that can be used for storing neuroimaging datasets (e.g. MRI, functional MRI, etc.). This format is adapted from the widely used ANALYZE 7.5 format from Mayo Clinic. A Nifti file consists of a header part containing the meta-data and the raw image data itself. The header and the image data can be stored together in the same file with an *.nii extension or can be split into two files with an *.hdr extension for the header and an *.img extension for the image data.

For functional MRI, it is not uncommon to store the dataset in a single 4D (3D space + time) Nifti file, especially when experiments involved several hundreds of image volumes. Just imagine the convenience of working with a single file as compared to several hundreds of files. However, there are also instances when it is convenient to manipulate datasets one volume at a time. In this case, having a single volume stored in a single file is very handy.

In this post, I will outline a simple Matlab script that can be used to extract the 3D images from the 4D NifTI file. The script uses some SPM functions such as spm_select(), spm_vol(), spm_read_vols(), and spm_write_vol(). Type “help function_name” in Matlab for a description of each function. Be sure that SPM is in Matla’s path before using this script. Without further ado, here is the full function list:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

function abk_4Dto3D(fname)
% function abk_4Dto3D(fname)
%
% Splits a 4D nifti file into a series of 3D nifti files.
% Needs SPM functions to work. If input file is fdata.nii,
% the output files will have filenames like fdata_001.nii,
% fdata_002.nii, etc.
 
if (nargin < 1)
    [fname,sts] = spm_select;
    if (sts == 0)
        fprintf('abk_4Dto3D: Operation cancelled.\n');
        return;
    end
end
 
vol = spm_vol(fname);
img = spm_read_vols(vol);
sz = size(img);
 
tvol = vol(1);
tvol = rmfield(tvol,'private');
tvol.descrip = 'generated by abk_4Dto3D.m';
 
[dn,fn,ext] = fileparts(fname);
 
for ctr=1:sz(4)
    tvol.fname = sprintf('%s%s%s_%.3d%s',dn,filesep,fn,ctr,ext);
    fprintf('Writing %s\n',tvol.fname);
    spm_write_vol(tvol,img(:,:,:,ctr));
end
fprintf('done.\n');
end

The script accepts a filename, given by the fname input parameter, of the 4D Nifti file. If fname is not specified, the script will prompt the user to select a file using spm_select(). Once the filename is specified, it will then read the 4D image data, and save it back one 3D image at a time.

Now to the code. Line 1 is simply the usual Matlab function declaration. This is followed by several lines of comments starting with the % character (lines 2 – 7). In line 9, the script checks if the function was called with an input parameter. If none was provided, it prompts the user to make a selection (line 10). It then checks if the user cancelled the operation (lines 11 – 14) and returns if true. Otherwise, it reads the header information using SPM’s spm_vol() function (line 17), loads the 4D data into img variable (line 18) using spm_read_vols(), and gets the image dimension (line 19).

In lines 21 – 23, the header information is copied into the tvol variable (line 21), which will be used later to save the images. The private field of tvol is then removed (line 22) and the header description is changed (line 23) to indicate the new file generator. In line 25, the filename is split into three variables dn, fn, and ext correspoding to the directory location, the base filename, and the extension, respectively.

In lines 27 – 31 is where the 4D image data is actually split into a series of 3D image files. Note how the variable tvol is repeatedly used, changing only the fname field each time. This is because the 3D images have basically the same header information.

Note that the same function could also be used to extract only certain volumes from the 4D nifti files with a slight modification. To implement this, we need to change the function declaration in line 1 into:

function abk_4Dto3D(fname,idx)

where the additional parameter idx is a vector containing the volume numbers to extract. We also need to change lines 27 – 31 as follows:

1
2
3
4
5

for ctr=1:length(idx)
    tvol.fname = sprintf('%s%s%s_%.3d%s',dn,filesep,fn,ctr,ext);
    fprintf('Writing %s\n',tvol.fname);
    spm_write_vol(tvol,img(:,:,:,idx(ctr)));
end

ctr now runs from 1 to the total number of volumes to extract given by the length of idx. The filenames will still be consecutively numbered from 1 to the length of idx (line 2). But this time, the 3D images that are saved are only the volumes specified in idx (line 4).

Have questions, leave a comment.

Re-posting this article on real-time functional MRI I wrote several years ago.

The technique known as functional magnetic resonance imaging (fMRI) has been extensively used to elucidate the functions of the human brain. Functional MRI provides a “window” where we can see what part of the brain is involved when we think, smell, taste, feel, or move. These windows are the activation maps indicating sites in the subject’s brain that are activated while the subject performs a given task and are usually obtained after analyzing voluminous amount of functional MRI images.

In a typical fMRI scanning session, the accumulated data can reach several hundreds of brain image volumes. The computational demand is usually very high that functional analyses are often delayed and activation maps are obtained long after the scanning session is completed and the subject is no longer inside the MR scanner. In most cases, this delay in analysis is often not an issue. However, the problem arises when after the analysis the data turned out to be corrupted with significant noise that no reliable activation could be extracted resulting in a waste of scanning time and delaying progress in research. To avoid this, the ideal scenario would be to “know” the results while the subject is still inside the MR scanner and before the experiment is concluded. This scenario calls for the immediate analysis of fMRI images as they are acquired and thus motivated the development of the real-time analysis of fMRI data.

Fig. 1. The conceptual framework of the real-time functional magnetic resonance imaging system composing of an MR scanner, a data storage device, and a computational server.

The real-time analysis of fMRI data offers several advantages both for theoretical studies and clinical applications. With real-time analysis, monitoring the task performance of the subject and the resulting quality of the acquired data can be easily achieved. It could also make functional mapping experiment more interactive by allowing ongoing paradigms to be altered if a need arises, making fMRI a more flexible tool for neurological investigations. Moreover, it also enables researchers to locate regions of interest that could be used for the next experimental run.

Fully real-time fMRI analysis could also provide an immediate feedback of the subject’s “ongoing” brain activity and thus enabling researchers to investigate the dynamical nature of the human brain. From a clinical perspective, some applications of real-time fMRI could be in pre-surgical planning. For instance, a surgeon operating a lesion could use functional brain studies to minimize the extent of the damage that could result in the operation. Real-time results of these functional studies are therefore critical. In spreading disorders such as Jacksonian seizures or migraine, the possibility of observing activation maps of such phenomena in real time could lead to a better understanding of the spreading mechanism of the disorder as well as to the development of therapeutic interventions to arrest the symptoms progression.

By providing an immediate feedback of the patient’s brain activity, real-time fMRI could also be used in assessing recovery treatment after a loss of, say for example motor function. These and several other potential applications are the compelling factors that motivate us to develop a fully real-time analysis system for fMRI time series. In simple terms, the research’s goal is to see the activation map unfolds in time as the subject performs the designated task, thereby enabling us to look how the brain works in real time.

Fig. 2. The above schematic shows the computational flowchart of the real-time analysis of fMRI time series, which starts immediately after an image is acquired and ends after the statistical parametric map is updated. The steps are repeated when a new image data is acquired.

In order to achieve this goal, we need to overcome two main difficulties, namely: 1) the need for a general real-time parametric analysis tool and 2) the computational requirement. The former requires analysis tools that we can use in real-time to be able to process fMRI data in the same rigor as that of offline analysis methods, while the latter demands computational prowess that can cope up with the needed computations. To overcome the first problem, we develop an algorithm to estimate the coefficients of general linear model (GLM), a versatile tool for parametric analysis, using an orthogonalization procedure. The algorithm offers several advantages including incremental estimation, minimal use of memory during the estimation process, fixed computational cost for each estimate update, among others, making it highly suitable for real-time applications. Using this as the basis for general real-time parametric analysis, a real-time analysis system is developed. To achieve results of comparable quality to that obtained by offline processing methods, realignment for motion correction and smoothing operations are incorporated into the computational pipeline. To meet the computational demand, we employ inexpensive and economically viable cluster of personal computers or PC cluster. The result is the real-time fMRI analysis system schematically shown in Fig. 1.

Fig. 3. Real-time activation maps at different time steps (n = 30, 50, 70, 90, 110, and 130) of a slice containing the primary motor cortex.

The system is composed of an MR scanner subsystem for data acquisition and paradigm control, a computational server (a PC cluster) for real-time fMRI data analysis, and a storage device for storing data. For offline operation, the MR scanner subsystem sends data to the storage device where the computational server can access for future analysis. For real-time operation, the MR scanner subsystem sends the data directly to the computational server for immediate processing. With this system, we achieved a fully real-time analysis of high spatial resolution whole-brain fMRI time series. The analysis includes realignment, smoothing, GLM estimation, and statistics computation (Fig. 2). All computations were performed immediately after the acquisition of each image volume and completed within TR set to 3 s. Real-time activation maps at different times for a slice containing the primiary motor area are shown in Fig. 3.

The real-time analysis system is designed such that the computational server does not have to be in the same site as the MR scanner. The idea is for imaging facilities that do not have a dedicated computational facility to be able to use remote computing resources on demand. This means that when needed the MR scanner subsystem could connect to a remote computational server via a high speed network and still be able to perform the required real-time operations. One way to realize this is to employ GRID technology. This will be the next step we will be taking, that is, the development of GRID-enabled real-time fMRI analysis system. The motivation is to make imaging facilities real-time capable by simply connecting to the GRID.

Just want to share this YouTube video from MicrosoftResearch. It’s just amazing. The idea is to augment the real environment with a virtual one and interact with it. Using four ceiling-mounted Kinects, the entire room is mapped into a virtual 3D environment. Virtual objects can be placed in this virtual environment and projected into the real one using handheld projectors. The interaction is quite impressive.

This is a brief guide on how to setup a backup SVN repository using svnsync. I will assume that subversion is already installed and the master repository is located in the machine called remote.com and in the directory /home/ruser/master_repos directory.

1. Create the mirror repository mirror_repos/ using svnadmin command

[user@mycomp] svnadmin create mirror_repos

2. Go to hooks/ subdirectory in the newly created repository

[user@mycomp] cd mirror_repos/hooks

You should see several *.tmpl files. One of these files is pre-revprop-change.tmpl. Copy it to pre-revprop-change, that is, without the .tmpl extension

[user@mycomp] cp pre-revprop-change.tmpl pre-revprop-change

Open this file using any text editor. Change it to the following:

1
2
3
4
5
6

USER="$2"
 
if [ "$USER" = "username" ]; then exit 0; fi
 
echo "Only username can change revprops" >&2
exit 1

Change “username” to the actual user name. Save the file and exit. Make sure it is also executable.

[user@mycomp] chmod 755 pre-revprop-change

3. Initialize the mirror repository using the svnsync init command

[user@mycomp] svnsync init DEST_URL SOURCE_URL
Copied properties for revision 0

In this example, DEST_URL is file:///home/user/mirror_repos and SOURCE_URL is svn+ssh://ruser@remote.com/home/ruser/master_repos

This initialize the mirror repository and it is now ready to be populated.

4. Start synchronizing. Use svnsync sync command to do it:

[user@mycomp] svnsync sync file:///home/user/mirror_repos
ruser@remote.com's password: xxxxxx
Committed revision 1.
Copied properties for revision 1.
Transmitting file data ........
Committed revision 2.
Copied properties for revision 2.
Transmitting file data ....
:

In case you want to commit back to the master repository, you can do it by issuing the command svn switch –relocate FROM_URL TO_URL before committing the changes. To do this, you also need to have the same UUID between the mirror and master repository.

Get the UUID from the master repository and copy it to the mirror repository. In the system where the master repository is running, run svnadmin dump

[ruser@remote] svnadmin dump -r0 /home/ruser/master_repos | head -n 3 > saved-uuid

Load it in the mirror repository

[user@mycomp] svnadmin load --force-uuid /home/user/mirror_repos < saved-uuid

The mirror and master repositories should now have the same UUID.

That’s it!

Definition

Given two random variables X and Y, Pearson’s correlation coefficient ρ is defined as the ratio between the covariance of the two variables and the product of their standard deviations:

where are the means and standard deviations of X and Y. The correlation coefficient is usually used as a measure of the strength of the linear relation between the two variables. Substituting the values of the covariance and standard deviations computed from sample time series gives the sample correlation coefficient, commonly denoted as r:

where

Alternatively, r can also be written as

Testing for significance

To test for the significance of the estimated correlation coefficient against the null hypothesis that the true correlation is equal to 0, one can compute the statistic

which has a Student’s t-distribution in the null case (zero correlation) with N-2 degrees of freedom. For instance in Matlab, you can compute the p-value using the function tcdf(). In particular, p = 2*tcdf(-abs(t), N-2) will give you the p-value for a two-tail t-test for a given t.

Alternatively, one can also convert the correlation coefficient using Fisher transform, given by

approximately follows a normal distribution with mean and standard deviation . With this, a z-score can now be defined as

Under the null hypothesis that r = r₀ and given the assumption that the sample pairs are independent and identically distributed, z follows a bivariate normal distribution. Thus an approximate p-value can be obtained from a normal probability table.

One can also use the Fisher transformation to test if two correlations r₁ and r₂ are significantly different by computing the z-score using the formula:

which is distributed approximately as N(0,1) when the null hypothesis is true. Here and are the number of samples used to compute and , respectively.

Incremental algorithm to compute the correlation coefficient

Define as follows:   and  The correlation coefficient can now be written as

To estimate the correlation coefficient incrementally, use the following algorithm:

Initialize n = 1:

for n = 2 to N, compute

Re-compute r using the above equation

end

Listed below are definitions of some vector or matrix products. I will be needing these definitions in the succeeding posts. Thanks to Jetpack‘s Latex support, I can now include equations in my post.

Outer Product

Given a row vector  with n elements and another vector  with m elements, the outer product of u and v is given by

The outer product as defined above can also be written as a matrix multiplication:

Kronecker Product

Given two matrices A which is m-by-n and B, which is q-by-p, the Kronecker product C of A and B is given by:

where C is an mp-by-nq matrix.

Khatri-Rao Product

Using the above definition, we can now define the Khatri-Rao product C of two matrices A and B as a column-wise Kronecker product of the two matrices. First, rewrite the matrices A and B as follows:

and

where represent the columns of A and B, respectively. The Khatri-Rao product is given by

Whew! I thought including equations was easy. Just hover your mouse over the equations above and you’ll see how complicated the code is to write such simple equations. There should be a better way.  Anyone?

There are several ways to do this. One is GUI-based. In Fedora Core, use the Places menu in the top panel. Within this menu item, click the ‘Connect to Server…’ sub-menu. This will popup the Connect to Server dialog box shown below.

In the Service type, select Windows share. Then enter the Windows server name in the Server text field. You can also fill in the other optional information, such as the Share, Folder, User Name, and Domain Name fields. After providing the needed information, click the Connect button. If the share requires a password, you will be prompted to enter one in a dialog box shown below:

Enter the password and select how long will you want the password to last, then click Connect. If successful, you will be connected to the specified Windows share. If not, an error will occur. Check if all the specified information are correct. Consult your system admin if you don’t know some of the required information.

The other alternative is to automount the share every time you login. To do this, you need root permission in your Linux system. First, create the directory where you want to mount the share (e.g., in /mnt/windows):

[root@mycomp] mkdir /mnt/windows

Then edit the /etc/fstab file and add the following lines.

//myserver/myshare   /mnt/windows   cifs   username=myname,password=mypass,uid=myuid,gid=mygid 0 0

Replace myname and mypass with the actual username and password for your Windows share, and myuid and mygid with the actual Linux user id (or username) and group id. You can also specify a full path to myserver, say, //myserver.mywindows.com/. Save fstab. Type

[root@mycomp] mount /mnt/windows

in your shell to immediately mount the share. If this works, the next time you login, the share will be automatically mounted.

That’s it. Have questions, drop a comment below.

I have been using Matlab for almost all of my computing needs. It is very easy to program especially when dealing with matrices. It has also a rich library of graphical interface. You can easily build a highly interactive graphical user interface using its GUI development environment or GUIDE. If you are familiar with event-based programming, you will have no difficulty using GUIDE.

Here is a quick tip you might find useful when programming in Matlab.

Dual-Monitor Setup

If you have a dual-monitor setup (using two monitors to extend the screen) and wanted to display all graphics-related output in the other monitor, you can use the ‘MonitorPositions’ property of the root window (0) to get the coordinates and dimensions of the two monitors as seen by Matlab. The code will look like this:

>> pos = get(0,'MonitorPositions');

If Matlab detects your dual-monitor setup, pos is a 2×4 matrix with the first row corresponding to the coordinates and dimension of your first monitor and the 2nd row that of the second monitor. The first two elements is the x and y position, while the last two correspond to the width and height. So for example, to position your figure in the second monitor in full screen, you can simply call the figure command with the position property set to pos(2,:), that is,

>> hfig = figure('Position',pos(2,:));

To make your application robust, you must first check if you have a second monitor before using the above code or else you will get an ‘Index exceeds matrix dimensions‘ error. You can accomplish this using the following code:

pos = get(0,'MonitorPositions');
sz = size(pos);
if (sz(1) > 1)
  hfig = figure('Position',pos(2,:));
else
  hfig = figure('Position',pos);
end

This will display your figure window full screen in the second monitor, if you have one. If you don’t, this will still display the figure, full screen, in the current monitor. And that’s it.

Been out for a while. My last entry was last Nov 2009! It’s more than a year and a half ago! It’s time for a reboot. There are so many things to learn, so many ideas to explore. Time for an upgrade!

First task is to update to the latest version of WordPress. As of this writing, the most recent is WordPress 3.2.1. Dashboard -> Updates -> Update automatically -> Done! WordPress’ automatic update is a snap. Next is the look and feel of the blog. Time for a new theme. After a short search, I found Nikolaj Masnikov’s Emerald Stretch, a flexible-width, three-column theme, cool to the eye and liked it. Clicked install -> Activate -> Done! And my blog turned green.

Several plugins also need some updates. Click plugins -> Select plugins to be updated -> Automatic update -> Done! Almost there. I just need a few more plugins to install. Facebook integration will be nice. Search -> Facebook-> Submit. Out of several FB-related plugins, I decided to use Otto’s Simple Facebook Connect plugin, including the Comments and Like Button components. How about site stats? There are several plugins for that too. For this, I chose Automattic’s Jetpack plugin, which also includes several other features! It enables WP.me short links and even supports LaTex. Now I can include in my post mathematical equations nobody would care to understand. Isn’t it smart.

That’s it, time to reboot!

Youtube Direct is a new tool from Google/Youtube that makes it easier for media organizations or websites to request users to submit video clips, review submitted clips, and possibly re-broadcast these clips, just like CNN’s iReport.  It allows websites to embed Youtube’s upload functionality with ease, includes moderation panel to approve or deny submitted clips, and provides a link back to the website when videos are viewed on Youtube.

Here is a summary of my attempt to install Youtube Direct (YTD). YTD is designed to run  on Google App Engine (GAE) so you need to have a GAE account. Don’t worry, it’s free. After getting one, I downloaded the source and prepared it to be deployed in GAE. To do this, there are several things that need to be configured. One is to have Eclipse and the related Google Plugin.

Installing Eclipse. Download Eclipse IDE for Java EE Developers from the Eclipse project website. The file is eclipse-jee-galileo-SR1-win32.zip, which is about 190MB. Extract the eclipse directory from the zip file. This directory contains the eclipse.exe executable, which runs Eclipse.

Install Google Plugin for Eclipse. Run Eclipse. Then select the Help menu and choose Install New Software.. menu item. In the Work with text box, enter the following URL: http://dl.google.com/eclipse/plugin/3.5. Click the Add button and then the OK button, keeping the name blank as it will be retrieved automatically from the site. Next, click the triangle next to “Plugin” and “SDKs”. Check the boxes next to “Google Plugin for Eclipse 3.5” and “Google App Engine Java SDK.” You can also include the “Google Web Toolkit SDK” if you like. Click the Next button and accept the terms of service to install the plugin. When the installation is complete, restart Eclipse and the plugin is installed.

Install subclipse to manage SVN checkout. Subclipse can be found here. To install it in Eclipse, open Eclipse and click the Help menu. Select Install New Software menu item. In the Work with text box, enter the following URL: http://subclipse.tigris.org/update_1.6.x. Click Add button. Follow the same process as in Step 2. This time selecting Subclipse from the available options. Restart Eclipse to complete the installation.

Download the source code. Now, you are ready to download the YTD source code. In Eclipse, select File menu, then click the Import menu item. In the Import dialog box, click the small triangle next to SVN, and then select Checkout Projects from SVN. Click Next button. From the Checkout from SVN dialog box, choose Create a new repository location. Click Next. Specify the URL as follows: http://youtube-direct.googlecode.com/svn/tags/20091113. Click Next. Select the folder to be checkout by highlighting http://youtube-direct.googlecode.com/svn/tags/20091113. Click Finish. A new project called ytd[tags/20091113] should appear in the Project Explorer after Eclipse completed downloading the source code.

Copy war/WEB-INF/appengine-web_DEFAULT.xml into war/WEB_INF/appengine-web.xml. Edit this file. Put within the <application> tag the registered App Engine instance name of your application. Specify the version of the application within the <version> tag. Replace the value of com.google.ytd.YTDeveloperKey property with your registered Youtube developer key. If you don’t have a developer key, get one here.

Run the application by clicking the Run action in Eclipse. If everything is working fine, you should have a server running at http://localhost:8080/. While still in Eclipse, deploy the application to GAE by clicking the Deploy App Engine Project button in the toolbar.

At first, I encountered several errors. I didn’t have a clue why. When I tried to deploy the application to GAE, it failed because it cannot find javac in the path. As it turned out, I needed to set Java in Windows –> Preferences –> Java –> Installed JREs to the installed Java SDK, not just the JRE. Any way, after setting Java to the SDK,  all the errors disappeared and the application worked and I was able to deploy it to GAE.

GAE dashboard. Go to the admin panel of YTD by visiting http://YOUR_APP.appspot.com/admin. From this panel, you can create new assignments, moderate submitted videos, etc. You can also get the codes that you can paste into a page in your website to enable video submission for each assignment.

That’s it! I still have to see how this works. You can also check the Youtube Direct Getting Started Guide for more details. Till next post.