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); |
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 |
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 |
There isn’t much magic here, just a simple idea that I use almost daily.
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.
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.
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…
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.
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 (↑)
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.









Recent Comments