eXTReMe Tracker
Dec 152014
 

Converting nifti files to dicom using nifti toolbox in matlab

clear; clc;
 addpath('/usr/local/matlab/toolbox/niftitools');
 subjectdir='/media/storage/SPECT_mapping/bbobby';
 cd(subjectdir);
 [dicomfile, dicompath] = uigetfile('*.dcm','Select DICOM header for metadata');
 metadata = dicominfo(fullfile(dicompath, dicomfile));
 metadataISAS=metadata; metadataMPR=metadata; metadataFUS=metadata;
 
 metadataMPR.SeriesDescription = 'SAG_MPRAGE_T1_Recon'; metadataMPR.SeriesNumber = 1;
 metadataISAS.SeriesDescription = 'ISAS Hyperperfusion'; metadataISAS.SeriesNumber = 2;
 metadataFUS.SeriesDescription = 'MPRAGE ISAS Fusion'; metadataFUS.SeriesNumber = 3;
 
 cd(subjectdir);
 isas = uigetfile('hyper*.nii','Select hyperpefusion file');
 structural = uigetfile('mprage*.nii','Select MPRAGE');
 
 isas = load_untouch_nii(isas); isas = isas.img; isas = int16(isas);
 structural = load_untouch_nii(structural); structural = structural.img;
 
 display('... Merging Hyperperfusion and MPRAGE');
 
 merged = structural;
 merged(isas>0) = max(max(max(structural)))+10;
 
 display('... Finished Merging')
 
 isasdir = fullfile(subjectdir,'isas');
 if exist(isasdir) ~= 7
     mkdir(isasdir);
     mkdir(fullfile(isasdir,'SPECT_ISAS')); spectdir = fullfile(isasdir,'SPECT_ISAS');
     mkdir(fullfile(isasdir,'structural')); structdir = fullfile(isasdir,'structural');
     mkdir(fullfile(isasdir,'fusion')); fusiondir = fullfile(isasdir,'fusion');
 end
 
 cd(isasdir);
 
 % dicomwrite(reshape(x4,[170,256,1,256]), 'SPECT_ISAS.dcm', metadata, 'CreateMode', 'copy');
 display('... Creating DICOMS');
 
 for slicenum = 1:size(isas,1)
     
     cd(spectdir); metadataISAS.InstanceNumber = slicenum;
     dicomwrite(fliplr(rot90(squeeze(isas(slicenum,:,:)))), sprintf('SPECT_ISAS_%03d.dcm',slicenum), metadataISAS, 'CreateMode', 'copy');
     
     cd(structdir); metadataMPR.InstanceNumber = slicenum;
     dicomwrite(fliplr(rot90(squeeze(structural(slicenum,:,:)))), sprintf('structural_%03d.dcm',slicenum), metadataMPR, 'CreateMode', 'copy');
     
     cd(fusiondir); metadataFUS.InstanceNumber = slicenum;
     dicomwrite(fliplr(rot90(squeeze(merged(slicenum,:,:)))), sprintf('FUSION_MPR_ISAS_%03d.dcm',slicenum), metadataFUS, 'CreateMode', 'copy');
     
     display(sprintf('... Creating slice %d',slicenum))
 end
 
 cd(subjectdir);
Oct 232014
 

To compute a transpose of an input file, i.e. converting columns to rows and rows to column using bash, you can accomplish that using awk and while loop.

#!/bin/bash
cols=`cat $1 | head -1 | wc -w`;
for i in $(seq 1 $cols);
do
awk -v "n=$i" '{print $n}' $1 | tr '\n' ' ';  # variable column output
printf "\n";
done

This example uses a variable to output a specific column from the input.

awk -v "n=$variable" '{print $n} ' input.txt

Will pass the variable to n and the nth column gets printed from input file.

Oct 242012
 

1. Convert the brain.mgz volume to analyze:

mri_convert $SUBJECTS_DIR/bert/mri/brain.mgz brain.img

2. Align brain.img with your functional data in spm (eg, func.img)

3. Create your spm statistical output (eg, t.img)

4. Create the registration file:

tkregister2 --s subjectname --mov meanswaf.img --regheader --reg  register.dat --surf orig

This will bring up the tkregister window with the orig volume. Hit the “Compare” button to see the functional. The green line will be the surface. Make sure the alignment is good (ie, the green line follows the bright intensity patterns in the functional). Hit the “Save” button to save the registration. This will create a file called register.dat. If you make modifications to the registration and then want to view or edit it later, re-run the above command WITHOUT –regheader. This requires a lot of manual tweaking to get the functional brain aligned with Bert brain.

5. View your functional data in the volume:

tkmedit subjectname orig.mgz -overlay spmT_0001.img -overlay-reg register.dat

6. View your functional data on the surface:

tksurfer subjectname lh inflated -overlay spmT_0001.img -overlay-reg register.dat

Adopted from spmPainting

Sep 192012
 
Simple image compression algorithm using plain FFT computation

FFT2 based image compression

Test code to generate a circular mask that removes the low frequency components from an image.

Step through radius to get the percentage compression.

IM=imread('images_01.jpg');IM=IM(:,:,2); %Select a channel
%
cx=size(IM,2)/2; cy=size(IM,1)/2;
ix=size(IM,2);iy=size(IM,1);r=833;
[x,y]=meshgrid(-(cx-1):(ix-cx),-(cy-1):(iy-cy));
c_mask=((x.^2+y.^2)>r^2);
%
filtIM = fft2(IM).*c_mask;
filtIM = ifft2(filtIM);
IM2=abs(double(IM)-abs(filtIM));
%
x=size(IM,1)*size(IM,2); y=size(find(c_mask==1)); y=y(1);
res=y/x*100;
%
figure(1); 
subplot(2,2,1); imshow(c_mask);
subplot(2,2,2); imagesc(IM2); axis image; axis off;
subplot(2,2,3); imshow(IM); title('Origianl image');
subplot(2,2,4); imagesc(abs(filtIM)); axis image; axis off;
title(sprintf('Percent original: %0.2f',res));
Nov 162011
 

On most Linux/Unix based systems, you can save your Matlab figure as a postscript instead of saving it as a png or jpg. One advantage of saving figures in this manner is that you can save multiple images as separate pages on that .ps document, whereas with png/jpg you have to save each image as a separate file.

To print the current image to postscript, all you need to run is:

print(‘-dpsc2′,’filename’);

To have multiple images appended within the postscript document, you need to use the -append switch.

print(‘-dpsc2′, ‘-append’, ‘filename’);

Your output postscript will have the name filename.ps

To convert the postscripts to pdf format, just run the ps2pdf command in the terminal window:

ps2fpdf filename.ps

The output of this would be filename.pdf

And you’re done.

Apr 142011
 

Assuming the nifti toolbox is in Matlab path, we can get the 91x109x91 mask to have the same dimensions as the normalized images generated with bounding boxes.

If we are making a mask for hippocampus, first we save that mask from WFU Pickatlas. Then to make it 79x95x68 voxels, run the following small script.

x=load_untouch_nii(‘hipp.img’);
x=x.img;xdim=[1:6 86:91]; ydim=[1:6 102:109]; zdim=[1:11 80:91];
origin=[40 57 26]; datatype=16;
x(xdim,:,:)=[]; x(:,ydim,:)=[]; x(:,:,zdim)=[];
nii=make_nii(x, [2 2 2], origin, datatype);
save_nii(nii, ‘boxedhippo.nii’)

We can then use these masks for signal extraction or any further processing.

Mar 302011
 

FSL: flirt is used to compute an intial affine normalization of the T1 weighted images; this is then fed to fnirt to compute the overall transformation. flirt is also used to register the EPI’s to the subject’s structural image. This was then used along with fnirt-s warp in applywarp to normalize the EPIs.

Script for normalizing the T1-weighted structurals to the template:

Registering T1-structural to MNI152
bet my_structural my_betted_structural
 flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm_brain -in my_betted_structural -omat my_affine_transf.mat
 fnirt --in=my_structural --aff=my_affine_transf.mat --cout=my_nonlinear_transf --config=T1_2_MNI152_2mm
 applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm --in=my_structural --warp=my_nonlinear_transf --out=my_warped_structural

Registering functional data to MNI152 (via structural scan)

bet my_structural my_betted_structural
 flirt -ref my_betted_structural -in my_functional -dof 7 -omat func2struct.mat
 flirt -ref ${FSLDIR}/data/standard/MNI152_T1_2mm_brain -in my_betted_structural -omat my_affine_transf.mat
 fnirt --in=my_structural --aff=my_affine_transf.mat --cout=my_nonlinear_transf --config=T1_2_MNI152_2mm
 applywarp --ref=${FSLDIR}/data/standard/MNI152_T1_2mm --in=my_functional --warp=my_nonlinear_transf --premat=func2struct.mat --out=my_warped_functional
Mar 142011
 
#!/bin/bash
 #$x would be the input filename and %3d will pad it to be a 3 digit number with zero pads, %4d will make it 4 digit and so on.
 num=`expr match "$x" '[^0-9]*\([0-9]\+\).*'`
 paddednum=`printf "%03d" $num`
 echo $paddednum
Dec 012010
 

Script to preprocess 4D fMRI time series:

spm_defaults;
 global defaults;
%% Slice Timing
display('Step 1: Slice Timing');
sliceorder = [1:2:17 2:2:16];
 refslice = 17;
 TR = 1;
 precision = 5;
 nslices = length(sliceorder);
P = dir('4Dout*.nii');
 P = char(P.name);
 TA = TR - TR/nslices;
if precision > 0
 TAstr = num2str(TA, precision);
 TA = str2double(TAstr);
 end
timing(1) = TA / (nslices -1);
 timing(2) = TR - TA;
 spm_slice_timing(P, sliceorder, refslice, timing);
%% Realignment and Reslicing
display('Step 2: Realignment');
 P=dir('a4Dout*.nii');
 P=char(P.name);
 FlagsC = struct('quality',1,'fwhm',5,'rtm',1,'interp',1);   %Flags for realignment
 spm_realign(P,FlagsC);
display('Step 3: Reslice');
 %Flags for reslicing
 FlagsR = struct('mask', 1, 'mean', 1,'interp',4,'which',0, 'wrap', defaults.realign.write.wrap);
 spm_reslice(P,FlagsR);
 display('... saving motion parameters as a motion_report.pdf ');
 saveas(gcf, 'motion_report', 'pdf');
 close all;
%% Normalization
display('Step 4: Normalization');
template = fullfile(spm('Dir'),'templates', 'EPI.nii');
 srcimage = 'mean*.nii';
meanimg = dir(srcimage);
 P = char(meanimg.name);
 snmat_name = strrep(P, '.nii', '_sn.mat');
 if exist(snmat_name, 'file')
 disp(['... removing existing file named: ' snmat_name]);
 delete(snmat_name);
 end
 objmask_name = '';
spm_normalise(template, P, snmat_name, defaults.normalise.estimate.weight, objmask_name,defaults.normalise.estimate);
%% Write Normalized Images
display('Step 5: Writing Normalized Images');
PP = dir('a4D*.nii');
 PP = char(PP.name);
 spm_write_sn(PP, snmat_name, defaults.normalise.write);
%% Smoothing
display('Step 6: Smoothing ... ');
PP = dir('wa4D*.nii');
 PP = char(PP.name);
 PP = cellstr(PP);
for list = 1:length(PP)
 Q = strcat('s', PP{list});
 spm_smooth(PP{list}, Q, defaults.smooth.fwhm);
 end
 %display(sprintf('\nPreprocessing compled.\n\nTime taken %d sec.\n', x));
Getting data extracted from time series:
Create ROI maps with WFU Pickatlas and MarsBar.
 Extract the contents of masks by:
 XYZ=struct(roi); XYZ=XYZ.XYZ;
%% Convert XYZ to mm
XYZ=[XYZ; ones(1, size(XYZ,2))]; XYZmm=roi.mat*XYZ; XYZmm(end,:)=[];XYZmm=XYZmm';
%%
p=spm_vol('swaf00xx.nii');
 [time_data]=spm_get_data(p, XYZ);
 time_data=mean(time_data,2);
Mar 082010
 

In order to convert files in .img and .hdr format to a single .nii file, we first need to install the nifti tools toolbox for matlab.

nifti tools

Once the path for nifti tools has been set in matlab, the files in .hdr and .img format can be combined to create a .nii file by running the following two commands:

nii=load_nii(‘filename_with_hdr_extension’);
save_nii(nii, ‘desired_filename_for_niftifile.nii’);

If you have freesurfer installed, you can convert between the files types with mri_convert by specifying the input (-it) and the output (-ot) type.

mri_convert -i name_of_input_file -it nifti1 -ot nii -o name_of_output_file

If you have functional dataset with several volumes over the course of your experiment, you can use the fslmerge function available in FSL to combine all the nifti files into a single 4D file.

fslmerge -t output_name.nii *.nii

Where *.nii would cover the entire set of nifti files if they are arranged sequentially.

Its not required to type in the extension .hdr or .img for the input image or output image if you have specified the input and the output types. By specifying the input type of nifti1, you tell the program that it will be a .hdr/.img pair, and specifying output type as nii goes for having your output with .nii extension.

To gzip the nifti file for use in freesurfer:

gzip filename.nii filename1.nii.gz

This will create a gzipped nifti file.

To automatically segment the structural volume:

recon_all -i filename1.nii.gz -subjid SUBJID -all

This will segment the structures containing in the structural file filename.nii.gz. Complete information on using recon-all for freesurfer can be found at FreeSurfer Wiki.

For some reason, Freesurfer seems to work well in tcsh. I’ve often encountered problems of weird nature when running it under bash. In order to run your autosegmentation with freesurfer type in tcsh at the terminal, since by default it is set to bash in most linux distros.

Mar 032009
 

I had to create a movie/animation out of a series of images that I had in my folder. There were several methods described on several forums which suggested using GIMP for creating GIF files from JPG or PNG files when one is using Linux. Now the problem was that I had 70 files in total and it would have taken me a long to time to add all those images to GIMP in order to create my animation file.

Then I read somewhere that we can actually create GIF images in command line mode using Imagemagick. But to do so you need to have Imagemagick installed on your system. To install it in Ubuntu, run:

sudo apt-get install imagemagick

This program has several features that would let you do all sorts of things with images. I was interested in converting a bunch of PNG files into an animated GIF file. To do so, we need to  rename the files in a way that they would lie in a sequence when arranged alphabetically (if you plan on doing it the easy way). Suppose you have files named slide_01.png, slide_02.png, slide_02.png….slide_xx.png, and you want to convert them to movie.gif, we run:

convert -delay 10 -loop 0 slide*.png movie.gif

The parameter delay inserts a desired delay between two consecutive slides. The number x used for delay inserts 10x milliseconds of delay between two frames. Loop parameter 0 makes it repeat infinitely.

If you had files with non-uniform names, then you need to input each of them in a sequence after the delay and loop parameters. Suppose you have file summer.jpg, fall, winter.jpg, fall.jpg and spring.jpg and you want to order them as fall, winter, spring and summer in the gif image seasons.gif with a 1 second delay between each of them, use:

convert -delay 100 -loop 0 fall.jpg winter.jpg spring.jpg summer.jpg seasons.gif

Here are some animations of numbers from 0 to 9 with varying delays.

The delay values specified the above cases were: 1, 5, 15, 25, 50 and 100.

gif animation with 10 millisecond delay

GIF animation with 10 millisecond frame delay

gif animation with 50 millisecond delay

GIF animation with 50 millisecond frame delay

gif animation with 150 millisecond delay

GIF animation with 150 millisecond frame delay

gif animation with 250 millisecond delay

GIF animation with 250 millisecond frame delay

gif animation with 500 millisecond delay

GIF animation with 500 millisecond frame delay

gif animation 1 second delay

GIF animation 1 second frame delay

=================

Feb 182009
 

This evening I took a snapshot of my desk (which holds my desktop) with my cellphone camera. After transferring it via bluetooth, I opened the image and it was rather amusing to see the exact image of desk inside of an identical image. So just out of curiousity I snapped another photo with roughly the same FOV as before and loaded it on the computer, and again took a picture of that.

Photo iterations: or photo within photo

Photo iterations: or photo within photo

The picture looks quite interesting, as it depicts the same object at different time points within the same spatial space. I don’t see any artistic value in this picture, but what I do see is a lot of potential for philosophical discussion about time and space and how we can represent variations in time within a 2D space.