Research: Math, Computing and MRI

Feb 082010

writeanalyze in Matlab

I am always look for different MRI file readers and writers for the myriad of formats that we use in MRI research. One of the relatively simple and common ones is the Analyze fileformat. Some of the large packages have writers (e.g., SPM) but I am typically wanting to do my own small processing and then write out the data. So, I wrote up my own writeanalyze.m function. It will do the basic formatting though the offsets etc don’t work. Try it out but I can’t guarantee anything.

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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104

function [] = writeanalyze(fname, data, ftype)

if( nargin == 2 )
ftype = 'int16';
end

if( strcmp( ftype, 'int16' ) == 1 )
file_type = 4; bpp = 16;
elseif( strcmp( ftype, 'uint16' ) == 1 )
file_type = 4; bpp = 16;
elseif( strcmp( ftype, 'int32' ) == 1 )
file_type = 8; bpp = 32;
elseif( strcmp( ftype, 'float' ) == 1 )
file_type = 16; bpp = 32;
elseif( strcmp( ftype, 'double' ) == 1 )
file_type = 64; bpp = 64;
else
error(sprintf('Unknown data type %s', ftype));
end

fp = fopen([fname '.hdr'], 'wb', 'b');

%%
%% Write the header_key part
%%
fwrite(fp, 348, 'int32');
fwrite(fp, repmat(' ', 1, 10), 'char');
fwrite(fp, repmat(' ', 1, 18), 'char');
fwrite(fp, 16384, 'int32');
fwrite(fp, 0, 'int16');
fwrite(fp, 'r ', 'char');

%%
%% Write the image_dimension part.
%%
fwrite(fp, length( size(data) ), 'int16');
for ii=1:length( size(data) )
fwrite(fp, size(data,ii), 'int16');
end

for ii=length( size(data) )+1:7
fwrite(fp, 1, 'int16');
end

fwrite(fp, 0, 'int16'); % unused 8
fwrite(fp, 0, 'int16'); % unused 9
fwrite(fp, 0, 'int16'); % unused 10
fwrite(fp, 0, 'int16'); % unused 11
fwrite(fp, 0, 'int16'); % unused 12
fwrite(fp, 0, 'int16'); % unused 13
fwrite(fp, 0, 'int16'); % unused 14

% data type
fwrite(fp, file_type, 'int16'); % 4 = signed short
fwrite(fp, bpp, 'int16'); % bpp
fwrite(fp, 0, 'int16');
for ii=1:8
fwrite(fp, 1.0, 'float32');
end
fwrite(fp, 0, 'float32');
fwrite(fp, 0, 'float32'); % funused 1
fwrite(fp, 0, 'float32'); % funused 2
fwrite(fp, 0, 'float32'); % funused 3

fwrite(fp, max(data(:)), 'float32');
fwrite(fp, min(data(:)), 'float32');
fwrite(fp, 0, 'float32');
fwrite(fp, 0, 'float32');

fwrite(fp, round(max(data(:))), 'int32'); % glmax
fwrite(fp, round(min(data(:))), 'int32'); % glmin

%%
%% Data history
%%
fwrite(fp, repmat(' ', 1, 80), 'char'); % descrip
fwrite(fp, repmat(' ', 1, 24), 'char'); % aux_file
fwrite(fp, '3', 'char'); % aux_file
fwrite(fp, repmat(' ', 1, 10), 'char'); % originator
fwrite(fp, repmat(' ', 1, 10), 'char'); % originator
fwrite(fp, repmat(' ', 1, 10), 'char'); % originator
fwrite(fp, repmat(' ', 1, 10), 'char'); % originator
fwrite(fp, repmat(' ', 1, 10), 'char'); % originator
fwrite(fp, repmat(' ', 1, 10), 'char'); % originator
fwrite(fp, repmat(' ', 1, 3), 'char'); % originator

fwrite(fp, 0, 'int32'); % views
fwrite(fp, 0, 'int32'); % vols_added
fwrite(fp, 0, 'int32'); % start_fiedl
fwrite(fp, 0, 'int32'); % field_skip
fwrite(fp, 0, 'int32'); % omax
fwrite(fp, 0, 'int32'); % omin
fwrite(fp, 0, 'int32'); % small_max
fwrite(fp, 0, 'int32'); % small_min

fclose(fp);

%%
%% Write the data
%%
fp = fopen([fname '.img'], 'wb', 'b');
fwrite(fp, data, ftype);
fclose(fp);

Feb 032010

A great comic

General, Math No Responses »

From xkcd:

Purity

Feb 012010

Reading Raw Data in Python

Computing, Python No Responses »

In a similar vein to reading raw data into Matlab, I created a similar type of function in Python:

1
2
3
4
5
6
7
8
9
10
11
12
13
14

def readraw(filename, shape, intype='int16', byteSwap=False):
""" readraw - To read in a raw file and reformat it to the right shape """

# Read in the file
if filename.endswith('gz'):
fp = gzip.open(filename, 'rb')
else:
fp = open(filename, 'rb')

d = fromfile(file=fp, dtype=intype).reshape(shape)

d.byteswap(byteSwap)

return d

Feb 012010

Reading raw data in Matlab

Computing, Matlab No Responses »

One of the most common things I do in Matlab almost always involves reading in binary data. For a few years I went through the typical fp=fopen('filename.dat'.... After typing the fopen, fread, reshape and fclose too many times, I finally made it into an all-in-one Matlab function called readraw() which will do all the reading and reformatting in one function call.

1
2
3
4
5
6
7
8
9
10
11

function [d] = readraw(filename, type, ds, endian)

if( nargin == 3 )
endian = 'b';
end

fp=fopen(filename, 'rb', endian);
d = fread(fp, prod(ds), type);
fclose(fp);

d = reshape(d, ds);

There isn’t much magic here, just a simple idea that I use almost daily.

Jan 272010

Huge Java Resource

Computing, Java No Responses »

I have been working on some numerical problems in the last few weeks. Mostly related to curve fitting and interpolation. I am slowly sliding back to Java bit-by-bit though I am not sure if I will give up Python as the syntax is so tight and it is actually very fast.

Anyway, I was looking around for some curve fitting and interpolation code on the Net and found Michael Flannigan’s website which has a great resource of math, stats, optimization and some other more subtle items. His Java code is here.

I can’t even do it justice by showing the list of packages as it is about twice as long as the screen capture shown to the right.

Jan 262010

Zotero – Reference Manager

General, Research No Responses »

There are many reference managers out there to use for organizing the myriad of references and papers one acquires during their career. One that I found recently is Zotero which is a javascript add-on for Firefox. I won’t go through all the interesting features here, but basically, it is extremely easy to save information you find on the Web. Not only web pages, but also journal article references (and PDFs) are easily saved. Another feature is one can sync the references to a database on their server and can be shared with groups of people.

Obviously a reference manager isn’t that useful if the information can’t be written into articles. So, the Zotero people have included mechanisms to pull references from their database and put it into documents being written in both Microsoft Word and OpenOffice. I have used it for one paper already and it worked out reasonably well.

So, how much does all this cost? Nothing, it is free. Try it out.

Jan 242010

How to organize a paper

General, Research No Responses »

When I first started in Science I can’t say anyone every sat down with me and told me how to write a paper. In reading papers, obviously something reasonably linear is easier. If one knows the background well, one can skip the Introduction and go right to the results. In writing a paper, I find it easier to think of it a little different. To that end I wrote up a template to help me order my thoughts and ideas. I will not say this it the only way to do it, nor the best. I find this way simple for me.

Below is the text. Attached here is a paper_outline one can print out.

The primary objective of this template is to help you organize your thoughts on the contents of a scientific paper. The standard format of a scientific paper is: Introduction, Methods, Results, Discussion and Conclusions. There are many modifications to this (e.g. inculde a Theory section or combine the Results and Discussion sections). Below is the likely order that one should work through to organize your paper, which can then be translated into the standard format (order). Try to move through each section in order, one follows naturally from the previous. Typically, each bullet point that you write below one of the sections is going to be a paragraph (though maybe not).

1 Primary Objectives

List the primary objective(s) of the work. Probably one primary objective and several secondary objectives.

2 Results to Fulfill Objectives

List the primary results that will fulfill the objectives. Will the results be shown by a figure, table or in the
text?

3 Discussion – Anomolies

List each discussion item. For each result, was there an anomoly about the result that needed to be described
further? Outliner in the data? Simple re-analysis with different starting conditions (e.g. curve fitting)?

4 Discussion – Pertaining to Previous Work

List each discussion item. For each result, how does it pertain to previous work? Are your results better? If
not, why?

5 Discussion – Future Work

List each extension to this work. How can the work be extended?

6 Describe the Methods

List each method that needs to be described. A method could be apparatus setup. Data analysis (and statistical tests). Participants (entry criterion, exclusion criterion). You will likely have several points here, again, each point will likely end up being a paragraph.

7 List the areas to Introduce

Given the results and discussion, what are the main areas that need to be introduced?

7.1 Motivation

7.2 Previous Work

7.3 Restate objectives

Jan 212010

Radiolab Podcast on Numbers

Math 2 Responses »

There is a great podcast (mp3 can be found on this page) from Radiolab that has a four part series on Numbers. The two parts that I enjoyed were on Benford’s Law (which I had not heard of before) and on the mathematician Paul Erdős.

Benford’s law basically states that the leading digit in real-life datasets is not uniformly distributed. Basically, the first digit of real-life datasets, the “1″ is more common than “2″ which is more common than “3″ etc. They give some background on how the law was found and interview someone using it in forensics.

Paul Erdős was a Hungarian born mathematician who is considered one of the brightest and most prolific mathematicians that has ever lived. The podcast gives some background to his life and the interactions he has with other mathematicians.

Very good podcast and well worth the listen…

http://www.wnyc.org/shows/radiolab/episodes/2009/10/09

Jan 192010

Anisotropic Diffusion Image Filtering in MRI

MRI, Matlab, Python No Responses »

Background

Magnetic resonance imaging has the tradeoff of signal-to-noise vs time vs resolution. You can only choose two. For some applications it may be better to get higher temporal and spatial resolution than signal-to-noise and then one may do some spatial filtering. Simple filtering would be applying a median filter or Gaussian smoothing over the image (or volume). But there are better techniques.

Smarter Filtering

One option for a smarter filter is the anisotropic diffusion filter which was first introduced to MRI in 1992 [1]. The basic idea is given a central voxel in a kernel and an estimation of noise the surrounding voxels are included in the smoothing based on the difference in signal to the central voxel relative to the estimation of noise.

I wrote a paper on this technique applied to multi-echo data [2].

There is a fine line between filtering and over-filtering. That is a whole separate discussion.

The images below are a single slice of an MPRAGE image without filtering (left) and with anisotropic diffusion filtering (right). The bottom set are just zoomed in versions of the top. The filtered data might be slightly over filtered but was done to show the affect of the filter.

Code

Matlab

The version below is for a 3D dataset:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

function [filt_vol] = aniso3d(orig_vol, kappa, niters)

if( nargin < 3 )
error('aniso3d: Need more parameters');
end

filt_vol = orig_vol;

for iters = 1:niters

dE = convn(filt_vol, [0 -1 1], 'full'); dE=dE(:,2:ncols(dE)-1,:);
dW = convn(filt_vol, [-1 1 0], 'full'); dW=dW(:,2:ncols(dW)-1,:);
dN = convn(filt_vol, [0; -1; 1], 'full'); dN=dN(2:nrows(dN)-1,:,:);
dS = convn(filt_vol, [-1; 1; 0], 'full'); dS=dS(2:nrows(dS)-1,:,:);
kernel = zeros(1,1,3); kernel(2) = -1; kernel(3) = 1;
dU = convn(filt_vol, kernel, 'full'); dU=dU(:,:,2:size(dU,3)-1);
kernel = zeros(1,1,3); kernel(1) = -1; kernel(2) = 1;
dD = convn(filt_vol, kernel, 'full'); dD=dD(:,:,2:size(dD,3)-1);

filt_vol = filt_vol + ...
3/28 * ((double(exp(- (abs(dE) / kappa).^2 )) .* double(dE)) - (double(exp(- (abs(dW) / kappa).^2 )) .* double(dW))) + ...
3/28 * ((double(exp(- (abs(dN) / kappa).^2 )) .* double(dN)) - (double(exp(- (abs(dS) / kappa).^2 )) .* double(dS))) + ...
1/28 * ((double(exp(- (abs(dU) / kappa).^2 )) .* double(dU)) - (double(exp(- (abs(dD) / kappa).^2 )) .* double(dD)));
end

For 4D data one can also smooth across the 4th dimension (whether it is time, diffusion etc).

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
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

function [filt_vol] = aniso3d_chan(orig_vol, kappa, niters)
%
% aniso3d_chan - Run the anisotropic diffusion filter in 3D
% and over the multiple channels.
%

if( nargin < 3 )
error('aniso3d: Need more parameters');
end

filt_vol = float(squeeze(orig_vol));

for iters = 1:niters
dE = convn(filt_vol, [0 -1 1], 'full'); dE=dE(:,2:ncols(dE)-1,:,:);
cE = repmat(sqrt(sum(dE.^2, 4)), [1 1 1 size(dE,4)]);
filt_vol = filt_vol + 3/28 * ((exp(- (cE / kappa).^2 )) .* (dE));
clear cE;
clear dE;

dW = convn(filt_vol, [-1 1 0], 'full'); dW=dW(:,2:ncols(dW)-1,:,:);
cW = repmat(sqrt(sum(dW.^2, 4)), [1 1 1 size(dW,4)]);
filt_vol = filt_vol - 3/28 * ((exp(- (cW / kappa).^2 )) .* (dW));
clear dW;
clear cW;

dN = convn(filt_vol, [0; -1; 1], 'full'); dN=dN(2:nrows(dN)-1,:,:,:);
cN = repmat(sqrt(sum(dN.^2, 4)), [1 1 1 size(dN,4)]);
filt_vol = filt_vol + 3/28 * ((exp(- (cN / kappa).^2 )) .* (dN));
clear dN;
clear cN;

dS = convn(filt_vol, [-1; 1; 0], 'full'); dS=dS(2:nrows(dS)-1,:,:,:);
cS = repmat(sqrt(sum(dS.^2, 4)), [1 1 1 size(dS,4)]);
filt_vol = filt_vol - 3/28 * ((exp(- (cS / kappa).^2 )) .* (dS));
clear cS;
clear dS;

kernel = zeros(1,1,3); kernel(2) = -1; kernel(3) = 1;
dU = convn(filt_vol, kernel, 'full'); dU=dU(:,:,2:size(dU,3)-1,:);
cU = repmat(sqrt(sum(dU.^2, 4)), [1 1 1 size(dS,4)]);
filt_vol = filt_vol + 1/28 * ((exp(- (cU / kappa).^2 )) .* (dU));
clear dU;
clear cU;

kernel = zeros(1,1,3); kernel(1) = -1; kernel(2) = 1;
dD = convn(filt_vol, kernel, 'full'); dD=dD(:,:,2:size(dD,3)-1,:);
cD = repmat(sqrt(sum(dD.^2, 4)), [1 1 1 size(dS,4)]);
filt_vol = filt_vol - 1/28 * ((exp(- (cD / kappa).^2 )) .* (dD));
clear dD;
clear cD;
end

Python

The Python code is very similar to the Matlab code above. It does 2D images or 3D volumes, but I have not coded the smoothing across the 4th dimension. That will have to be done later.

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

def aniso(v, kappa=-1, N=1):

if kappa == -1:
kappa = prctile(v, 40)

vf = v.copy()

for ii in range(N):
dE = -vf + roll(vf,-1,0)
dW = vf - roll(vf,1,0)

dN = -vf + roll(vf,-1,1)
dS = vf - roll(vf,1,1)

if len(v.shape) > 2:
dU = -vf + roll(vf,-1,2)
dD = vf - roll(vf,1,2)

vf = vf + \
3./28. * ((exp(- (abs(dE) / kappa)**2 ) * dE) - (exp(- (abs(dW) / kappa)**2 ) * dW)) + \
3./28. * ((exp(- (abs(dN) / kappa)**2 ) * dN) - (exp(- (abs(dS) / kappa)**2 ) * dS))
if len(v.shape) > 2:
vf += 1./28. * ((exp(- (abs(dU) / kappa)**2 ) * dU) - (exp(- (abs(dD) / kappa)**2 ) * dD))

return vf

References

1. G. Gerig et al., “Nonlinear anisotropic filtering of MRI data,” Medical Imaging, IEEE Transactions on 11, no. 2 (1992): 221-232. (↑)
2. Craig K Jones, Kenneth P Whittall, and Alex L MacKay, “Robust myelin water quantification: averaging vs. spatial filtering,” Magnetic Resonance in Medicine: Official Journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine 50, no. 1 (July 2003): 206-209 (↑)

Jan 182010

Pulse Sequence Diagrammer

MRI, Matlab No Responses »

Overview

The matlab code in this directory should facilitate creating publication quality PSDs (pulse sequence diagrams) using Matlab. Look at the example files (cse.m, cpmg.m and fse.m) to see how to use the code. All files are script files so this should run on any machine that Matlab runs on.

Code: mrpsd_12.tar.gz

Matlab Version …

I know that it worked under version 5.x of Matlab, but it should work under any newer version as well.

Why under Matlab?

Ahh.. good question. There are many reasons:

1) Many people use Matlab for their data analysis and general coding.

2) All of the print facilities are built in (so you can print to JPEG, Postscript, BMP, TIFF etc etc).

3) Many things come free with the way that it is designed, for example, if you want to look at only one temporal section of your PSD, all you have to do is plot it up and then do: set(gca, ‘xlim’, [50 100]) (if you want to look at between 50ms and 100ms). USE YOUR IMAGINATION HERE. There are potentially lots of little things like the previous example that I have not even thought of.

E-mail me

I would be very interested in any suggestions, fixes (!) that you can send along to make this toolbox better. I would also like any more example files that plot up other pulse sequences (spectroscopy, EPI etc etc). My e-mail is craig@mri.jhu.edu. It is free software, I will not restrict use in any way, shape or form (other than don’t sell it). I would appreciate, though, any enhancements that you can. I will try to make available updates as often as possible.

Standard Disclaimer

By using the software, I accept absolute no responsibility for anything. Use it at your own risk. It is absolutely GPL‘ed software.

Have fun with it.