# Huge Java Resource

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.

# Zotero – Reference Manager

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.

# How to organize a paper

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

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…

# Anisotropic Diffusion Image Filtering in MRI

## 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 ((G. Gerig et al., “Nonlinear anisotropic filtering of MRI data,” Medical Imaging, IEEE Transactions on 11, no. 2 (1992): 221-232. )).  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 ((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)).

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

## Code

### Matlab

The version below is for a 3D dataset:
[cc lang="matlab"]

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
[/cc]

For 4D data one can also smooth across the 4th dimension (whether it is time, diffusion etc).
[cc lang="matlab"]
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
[/cc]

### 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.
[cc lang="python"]
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
[/cc]

# Pulse Sequence Diagrammer

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

# Tracking Research

It can be an onerous task to keep track of all the research that gets published these days. I can’t say that I am great at it either, but there are some reasonably standard Google tools that I have found useful.  Back when I first started in MRI research, my supervisor had a lab meeting in which one person’s job (rotating, thankfully) was to go to the library and photocopy the Table of Contents of all the relevant journals, go through looking for relevant or interesting articles and present to the group.  Those were the days.

Google Reader is an RSS feed aggregator which basically means you can subscribe to RSS feeds and it gives an interface to be able to see them. Sites like PubMed are great for searching for journal articles (PubMed for the medical arena). Their site has changed a lot over the past 10 years. Recently, on the search results page there is an RSS feed symbol (orange symbol in the graphic above) that one can click. In most browsers this will bring up a window in which you can select how to watch the feed. There are lots of offline readers, but I like Google Reader.  Many journals also have RSS feeds of the articles that come out.

Another Google tool that works well is Google Docs. It is a good place to writeup documents that need to be shared with others and it keeps revisions of the documents. I can’t say I have used it for writing a paper, yet, but it is getting very close. The one thing that I haven’t worked out is how to put references in a Google Doc though I am sure there is a way.

# Noise in MRI (magnitude) data

## Background

Magnitude MRI data has Rician noise distribution by definition ((Hákon Gudbjartsson and Samuel Patz, “The Rician Distribution of Noisy MRI Data,” Magnetic resonance in medicine : official journal of the Society of Magnetic Resonance in Medicine / Society of Magnetic Resonance in Medicine 34, no. 6 (December 1995): 910-914)). It comes about because two channels each with Gaussian noise are squared and added together ((R M Henkelman, “Measurement of signal intensities in the presence of noise in MR images,” Medical Physics 12, no. 2 (April 1985): 232-233.)).  There is a longer description here.

## Modeling

The Rician noise is created as $latex y_e(t_i) = \sqrt{ \left[y(t_i) + e_1 \right]^2 + e_2^2 }$, where $latex y$ is the true signal, and $latex e_1$ and $latex e_2$ are random numbers from a Gaussian distribution with zero mean and standard deviation $latex \sigma$.  The standard deviation, $latex \sigma$, for the Gaussian distribution is related to the signal to noise ratio and is typically on the order of 1% – 10% of the signal $latex y$.

## Code

It is relatively easy to model this using Matlab or Python. For the code here I am modeling a T2 decay curve and then the noise.

### Matlab

[cc lang="matlab"]

%  Setup the initial variables
rho = 100;
t2 = 80; % in ms
te = 10:10:320;  % in ms

%  Create a T2 decay curve
y = rho * exp(-te / t2 );

%  Define the noise to be 5% of the signal
s = 5;

%  Create the two Gaussian random variable vectors
e1 = s * randn(size(y));
e2 = s * randn(size(y));

%  Now create the new, noisy decay curve.
y_e = sqrt( (y+e1).^2 + (e2).^2 );

[/cc]

### Python

The Python version is quite similar.

[cc lang="python"]

from __future__ import division

#  Setup the initial variables
rho = 100
t2 = 80 # in ms
te = r_[10:330:10] # in ms

#  Create a T2 decay curve
y = rho * exp( -te / t2 )

#  Define the noise to be 5% of the signal
s = 5;

#  Create the two Gaussian random variable vectors
e1 = normal(0, 5, y.shape)
e2 = normal(0, 5, y.shape)

#  Now create the new, noisy decay curve.
y_e = sqrt( (y+e1)**2 + (e2)**2 );

[/cc]

There are a couple of small gotcha’s that at least tripped me up as I am still relatively new to Python.

1. The first is that under Python 2.x all data is processed as integer (not doubles, as the default is in Matlab).  Supposedly this is going to change in Python 3, but to get around it for now, the best thing to do is to add the [cci lang="python"]from __future__ import divison[/cci].
2. To define [cci lang="python"]te[/cci] I had to go to 330, rather than 320 as the generator is an open set on the higher end so it does not include the number.
3. There are several options for creating the random numbers.  There is a Python module called [cci lang="python"]random[/cci] that could be used.  Instead I used the Numpy [cci lang="python"]normal[/cci] instead as I can pass in the shape parameter.

# everydns

Hmm… It seems that everydns.net has been sold to the dyndns people. I have used everydns.net for quite a few years now (even donated!). There are two things that I am nervous about doing myself, one is DNS and the second is MX (mail). I have been very happily pointing my MX records to Google and very happily using everydns.net to point my A and CNAME records where I want. I started looking through for a replacement (just in case, at some point they will likely start charging for the service). I found a few that seem free but then look like they are slightly limiting. One I found seems to be very similar to everydns.net (free and not limiting). One is called http://xname.org and another is http://www.freedns.ws.

# Change font size of tick labels

For making figures it is sometimes important (or quite important) to increase the font size of the x or y ticklabels. Here is one way I found to do it:

[cc lang="python"]
fig1 = figure()
for t in gca().get_yticklabels():
t.set_fontsize(14)

fig1.canvas.draw()

[/cc]

For some reason there has to be a [cci lang="python"]fig1.canvas.draw()[/cci] at the end of this to refresh the figure.