August 06, 2010, at 09:26 PM by pcook -

Changed line 11 from:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the available formats. The pointset2scheme? program converts point sets to scheme files. The point set for this data is stored in pigPoints.txt. The appropriate scheme file format depends on the application and the available information about the image acquisition. In this case, we know most of the imaging parameters, so we'll generate a detailed scheme file:

to:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the available formats. The pointset2scheme program converts point sets to scheme files. The point set for this data is stored in pigPoints.txt. The appropriate scheme file format depends on the application and the available information about the image acquisition. In this case, we know most of the imaging parameters, so we'll generate a detailed scheme file:

August 06, 2010, at 09:24 PM by pcook -

Changed line 34 from:

[@pdview -inputfile pig_dt_eig_sys.Bdouble -scalarfile pig_fa.hdr]

to:

pdview -inputfile pig_dt_eig_sys.Bdouble -scalarfile pig_fa.hdr

August 06, 2010, at 09:24 PM by pcook -

Changed lines 11-12 from:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the available formats. The pointset2scheme program converts point sets to scheme files. The point set for this data is stored in pigPoints.txt. The appropriate scheme file format depends on the application and the available information about the image acquisition. In this case, we know most of the imaging parameters, so we'll generate a detailed scheme file:

to:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the available formats. The pointset2scheme? program converts point sets to scheme files. The point set for this data is stored in pigPoints.txt. The appropriate scheme file format depends on the application and the available information about the image acquisition. In this case, we know most of the imaging parameters, so we'll generate a detailed scheme file:

Changed lines 32-36 from:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/fa05.png

to:

This command produces several Analyze images, and some Camino raw data files. All output is prepended with pig_. See the man page for analyzedti for more details. We can check the reconstruction using pdview:

[@pdview -inputfile pig_dt_eig_sys.Bdouble -scalarfile pig_fa.hdr]

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/pig_pdview.png

Changed line 42 from:

fa < PigBrainDT.Bdouble | shredder $((128*128*8*4)) $((128*128*8)) $((128*128*8*10)) > PigBrainFA_SLICE05.Bdouble \\

to:

fa < pig_dt.Bdouble | shredder $((128*128*8*4)) $((128*128*8)) $((128*128*8*10)) > PigBrainFA_SLICE05.Bdouble \\

Changed line 46 from:

shredder $((128*128*8*8*4)) $((128*128*8*8)) $((128*128*8*8*10)) < PigBrainDT.Bdouble | sfplot -xsize 128 -ysize 128 -inputmodel dt -backdrop PigBrainFA_SLICE05.Bdouble -dircolcode -maxd 0.4E-9 -projection 2 1 > PigBrain05_PAS.rgb \\

to:

shredder $((128*128*8*8*4)) $((128*128*8*8)) $((128*128*8*8*10)) < pig_dt.Bdouble | sfplot -xsize 128 -ysize 128 -inputmodel dt -backdrop PigBrainFA_SLICE05.Bdouble -dircolcode -maxd 0.4E-9 -projection 2 1 > PigBrain05_PAS.rgb \\

August 06, 2010, at 09:11 PM by pcook -

Changed lines 28-29 from:

ls *.hdr > imagelist.txt

to:

ls NEX*.hdr > imagelist.txt

Changed lines 35-70 from:

 
Alternatively, you can convert the camino data files to analyze format, just by adding a header, and view them in, for example, MRICro: 
 
analyzeheader -datadims 128 128 10 -voxeldims 0.5 0.5 0.5 -datatype double > fa.hdr \\

cp PigBrainFA.Bdouble fa.img

dteig computes the eigensystem of each tensor:

dteig < PigBrainDT.Bdouble > PigBrainEIG.Bdouble

We can generate colour-coded direction maps easily from the output in Matlab:

>> fid = fopen('PigBrainEIG.Bdouble', 'r', 'b');
>> data = fread(fid, 'double'); \\ >> fclose(fid);
>> % Reshape into slices; 12 values per voxel
>> eigs = reshape(data, 12, 128,128,10);
>> % Elements 2, 3 and 4 are the principal eigenvector
>> fapds = abs(permute(eigs(2:4,:,:,:), [2,3,4,1]));
>> % Weight with fa.
>> fapds(:,:,:,1) = fapds(:,:,:,1).*fa;
>> fapds(:,:,:,2) = fapds(:,:,:,2).*fa;
>> fapds(:,:,:,3) = fapds(:,:,:,3).*fa;
>> fapds=permute(fapds, [1,2,4,3]);
>> % Transpose to standard orientation.
>> for i=1:3;
>> for j=1:10
>> fapds(:,:,i,j)=flipud(fapds(:,:,i,j)');
>> end;
>> end
>> imshow(fapds(:,:,:,5))
>> for i=1:10
>> imwrite(fapds(:,:,:,i), sprintf('pds%02i.png', i));
>> end

That outputs the following image (slice 5 only):
http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/pds05.png

to:

Diffusion tensor plot

The command sfplot creates plots of diffusion MRI reconstructions over slices. In particular, it can create plots of the diffusion tensor over a slice. It can also create plots of other reconstructed features, such as the PAS or q-ball ODF; see for example the Parallel PASMRI case study or the sfplot man page. The following commands create a plot of the DT over slice 5 of the pig brain data set. We start by creating a fractional anisotropy map to use as background for the plot: \\

August 06, 2010, at 09:03 PM by pcook -

Changed lines 7-8 from:

The data set in this case study comes from a post-mortem porcine brain and was kindly provided by Tim Dyrby from The Danish Research Centre for Magnetic Resonance, Copenhagen University Hospital, Hvidovre, Denmark. The data set is intended for demonstration and testing of the Camino software package only. If you wish to use this data for other purposes, such as publications, other web sites or commercial enterprises, please contact Tim at (timd|at|drcmr.dk) before doing so. You can download the data set . The acquisition uses a spherical acquisition scheme with 61 unique gradient directions. The gradient directions come from the file camino/PointSets/Elec061.txt and minimize the electrostatic energy for 61 pairs of equal-and-opposite point electrical charges free to move on the surface of the sphere. The procedure outlined in [Jansons and Alexander, Inverse Problems, Vol 19, pp 1031-1046, 2003] computed them. The diffusion-weighting gradient strength G = 61 mT m^-1; pulse onset separation DELTA = 30ms; pulse width delta = 23ms; TE = 60ms; b = 3146 s mm^-2; NEX=2. The protocol also includes 3 acquisitions with no diffusion weighting. The image has 10 slices with in-plane resolution 128x128; voxel size 0.5x0.5x0.5 mm³. You can download the data here.

to:

The data set in this case study comes from a post-mortem porcine brain and was kindly provided by Tim Dyrby from The Danish Research Centre for Magnetic Resonance, Copenhagen University Hospital, Hvidovre, Denmark. The data set is intended for demonstration and testing of the Camino software package only. If you wish to use this data for other purposes, such as publications, other web sites or commercial enterprises, please contact Tim at (timd|at|drcmr.dk) before doing so. You can download the data set . The acquisition uses a spherical acquisition scheme with 61 unique gradient directions. The gradient directions come from the file camino/PointSets/Elec061.txt and minimize the electrostatic energy for 61 pairs of equal-and-opposite point electrical charges free to move on the surface of the sphere. The procedure outlined in [Jansons and Alexander, Inverse Problems, Vol 19, pp 1031-1046, 2003] computed them. The diffusion-weighting gradient strength G = 61 mT m^-1; pulse onset separation DELTA = 30ms; pulse width delta = 23ms; TE = 60ms; b = 3146 s mm^-2; NEX=2. The protocol also includes 3 acquisitions with no diffusion weighting. The image has 10 slices with in-plane resolution 128x128; voxel size 0.5x0.5x0.5 mm³. You can download the data here.

Changed lines 11-24 from:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the format. The simple perl script camino/PointSets/PointSetToScheme converts pigPoints.txt to a scheme file:

PointSets/PointSetToScheme pigPoints.txt 375330 3 0.0223

The first argument is the point set file; the second is the wavenumber q = gamma*delta*G = 2.6751987E8 * 0.023 * 0.061, where gamma is the gyromagnetic ratio of protons in water; the third argument is the number of b=0 images and the last is the diffusion time. Here we use the standard approximation of the diffusion time for finite pulse widths t = DELTA - delta/3 = 0.03 - 0.23/3. This is the resulting . A scheme file specifies the gradient directions along with the zero directions. The zero directions should always be positioned before the non-zero directions.

Data File

Each acquisition from the scanner is in an analyze format file NEX2_s_20060127_001_semsdw_15_image??, where ?? goes from 01 to 64; 01-03 are the b=0 images and 04-64 are the diffusion-weighted images with gradients in the same order as the scheme file. Each file stores image intensities as little-endian 2-byte short integers. We must generate a data file in camino format, which requires switching to big endian representation and combining all the analyze files into a single file with voxel ordering (see camino/man/man1/camino.1). Here is the command:

cat *.img | shredder 0 -2 0 | scanner2voxel -voxels 163840 -components 64 -inputdatatype short > PigBrain.Bfloat

The shredder command switches the byte ordering and scanner2voxel ensures that all measurements for each voxel are contiguous in the data file. The datafile PigBrain.Bfloat contains big-endian 4-byte floating-point data.

to:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the available formats. The pointset2scheme program converts point sets to scheme files. The point set for this data is stored in pigPoints.txt. The appropriate scheme file format depends on the application and the available information about the image acquisition. In this case, we know most of the imaging parameters, so we'll generate a detailed scheme file:

pointset2scheme -inputfile pigPoints.txt -version STEJSKALTANNER -dwparams 0.061 0.03 0.023 0.06 > PigBrain_ST.scheme

The "version" argument tells the program that we want to describe a Stejskal-Tanner sequence, which uses two rectangular gradient pulses of equal duration. The "dwparams" argument is followed by the parameters of the sequence: gradient strength (0.061 T / m), pulse separation (0.03 s), pulse duration (0.023 s), and echo time (0.06 s). Alternatively, we could compute the b-value for the scheme, 3146 s/mm^2, and created a less detailed scheme file:

pointset2scheme -inputfile pigPoints.txt -version BVECTOR -dwparams 3146 > PigBrain_BVEC.scheme

For most applications, different scheme files are equivalent as long as they describe the same b-value. Some programs, like the monte-carlo diffusion simulation, require a more detailed description of the pulse sequence.

Changed lines 24-68 from:

To fit the diffusion tensor in each voxel we can use dtfit or modelfit:

dtfit PigBrain.Bfloat PigBrain.scheme > PigBrainDT.Bdouble

or equivalently:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 > PigBrainDT.Bdouble

Both commands give the same output. Note that the default output data type is big-endian doubles. The modelfit program has more options (see the modelfit man page) than dtfit, such as thresholding the b=0 measurements to remove the noisy background you can see in many of the images below, like this:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 -bgthresh 500 > PigBrainDT.Bdouble

Simple tensor statistics

From the diffusion tensor, we can compute simple statistics such as Trace(D), fractional anisotropy and principal directions. The commands trd and fa compute the first two:

trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble
fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble

A simple way to look at the results is to load the data files into matlab and display with imshow (use image instead if you do not have the image-processing toolbox) or write out slice by slice:

>> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b');
>> data = fread(fid, 'double');
>> fclose(fid);
>> % Normalize and reshape into slices
>> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data));
>> % Show one slice
>> imshow(flipud(trd(:,:,5)'));
>> % Write out all slices as .png images
>> for i=1:10
>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));
>> end

>> % Repeat for FA
>> fid = fopen('PigBrainFA.Bdouble', 'r', 'b');
>> data = fread(fid, 'double');
>> fclose(fid);
>> fa = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data));
>> imshow(flipud(fa(:,:,5)'));
>> for i=1:10
>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));
>> end

That outputs the following images (slice 5 only): \\

to:

Each acquisition from the scanner is in an analyze format file NEX2_s_20060127_001_semsdw_15_image??, where ?? goes from 01 to 64; 01-03 are the b=0 images and 04-64 are the diffusion-weighted images with gradients in the same order as the scheme file.

To fit the diffusion tensor in each voxel we can use analyzedti:

ls *.hdr > imagelist.txt

analyzedti imagelist.txt ./pig_ PigBrain_ST.scheme

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/fa05.png \\

December 09, 2009, at 03:11 PM by shahrum - added title, link to pig brain data

Changed lines 1-2 from:

Tutorial: detecting crossing fibres

to:

(:title Tutorial: detecting crossing fibres:)

Changed line 7 from:

The data set in this case study comes from a post-mortem porcine brain and was kindly provided by Tim Dyrby from The Danish Research Centre for Magnetic Resonance, Copenhagen University Hospital, Hvidovre, Denmark. The data set is intended for demonstration and testing of the Camino software package only. If you wish to use this data for other purposes, such as publications, other web sites or commercial enterprises, please contact Tim at (timd|at|drcmr.dk) before doing so. You can download the data set . The acquisition uses a spherical acquisition scheme with 61 unique gradient directions. The gradient directions come from the file camino/PointSets/Elec061.txt and minimize the electrostatic energy for 61 pairs of equal-and-opposite point electrical charges free to move on the surface of the sphere. The procedure outlined in [Jansons and Alexander, Inverse Problems, Vol 19, pp 1031-1046, 2003] computed them. The diffusion-weighting gradient strength G = 61 mT m^-1; pulse onset separation DELTA = 30ms; pulse width delta = 23ms; TE = 60ms; b = 3146 s mm^-2; NEX=2. The protocol also includes 3 acquisitions with no diffusion weighting. The image has 10 slices with in-plane resolution 128x128; voxel size 0.5x0.5x0.5 mm³.

to:

The data set in this case study comes from a post-mortem porcine brain and was kindly provided by Tim Dyrby from The Danish Research Centre for Magnetic Resonance, Copenhagen University Hospital, Hvidovre, Denmark. The data set is intended for demonstration and testing of the Camino software package only. If you wish to use this data for other purposes, such as publications, other web sites or commercial enterprises, please contact Tim at (timd|at|drcmr.dk) before doing so. You can download the data set . The acquisition uses a spherical acquisition scheme with 61 unique gradient directions. The gradient directions come from the file camino/PointSets/Elec061.txt and minimize the electrostatic energy for 61 pairs of equal-and-opposite point electrical charges free to move on the surface of the sphere. The procedure outlined in [Jansons and Alexander, Inverse Problems, Vol 19, pp 1031-1046, 2003] computed them. The diffusion-weighting gradient strength G = 61 mT m^-1; pulse onset separation DELTA = 30ms; pulse width delta = 23ms; TE = 60ms; b = 3146 s mm^-2; NEX=2. The protocol also includes 3 acquisitions with no diffusion weighting. The image has 10 slices with in-plane resolution 128x128; voxel size 0.5x0.5x0.5 mm³. You can download the data here.

October 22, 2009, at 02:33 AM by shahrum -

Changed lines 151-167 from:

>> fid = fopen('PigBrainVC.Bint', 'r', 'b'); \\ >> vc=fread(fid, 'int'); \\ >> fclose(fid); \\ >> vc=reshape(vc, 128,128,10); \\ >> % Repeat for FA \\ >> fid = fopen('PigBrainFA.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> fa = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> for i=1:10 \\ im(:,:,1) = fa(:,:,i).*(1-0.5*(vc(:,:,i)>=4)) + 0.5*(vc(:,:,i)>=4); \\ im(:,:,2) = fa(:,:,i).*(1-0.5*(vc(:,:,i)>=4)) + 0.5*(vc(:,:,i)>=4); \\ im(:,:,3) = fa(:,:,i); \\ imwrite(im, sprintf('/tmp/favc%02i.png', i)); \\ end \\\\ Here is the image for slice 5: \\

to:

>> fid = fopen('PigBrainVC.Bint', 'r', 'b');
>> vc=fread(fid, 'int');
>> fclose(fid);
>> vc=reshape(vc, 128,128,10);
>> % Repeat for FA
>> fid = fopen('PigBrainFA.Bdouble', 'r', 'b');
>> data = fread(fid, 'double');
>> fclose(fid);
>> fa = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data));
>> for i=1:10
>> im(:,:,1) = fa(:,:,i).*(1-0.5*(vc(:,:,i)>=4)) + 0.5*(vc(:,:,i)>=4);
>> im(:,:,2) = fa(:,:,i).*(1-0.5*(vc(:,:,i)>=4)) + 0.5*(vc(:,:,i)>=4);
>> im(:,:,3) = fa(:,:,i);
>> imwrite(im, sprintf('/tmp/favc%02i.png', i));
>> end
Here is the image for slice 5: \\

October 22, 2009, at 02:30 AM by shahrum -

Changed line 127 from:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/PigBrain05_PAS.png \\

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/PigBrain05_PAS.png \\

October 22, 2009, at 02:28 AM by shahrum -

Changed line 142 from:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/VC_ThreshSelect.png \\

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/VC_ThreshSelect.png\\

October 22, 2009, at 02:13 AM by shahrum -

Changed line 19 from:

Each acquisition from the scanner is in an analyze format file NEX2_s_20060127_001_semsdw_15_image??, where ?? goes from 01 to 64; 01-03 are the b=0 images and 04-64 are the diffusion-weighted images with gradients in the same order as the scheme file. Each file stores image intensities as little-endian 2-byte short integers. We must generate a data file in camino format, which requires switching to big endian representation and combining all the analyze files into a single file with voxel ordering (see camino/man/man1/camino.1). Here is the command: \\

to:

Each acquisition from the scanner is in an analyze format file NEX2_s_20060127_001_semsdw_15_image??, where ?? goes from 01 to 64; 01-03 are the b=0 images and 04-64 are the diffusion-weighted images with gradients in the same order as the scheme file. Each file stores image intensities as little-endian 2-byte short integers. We must generate a data file in camino format, which requires switching to big endian representation and combining all the analyze files into a single file with voxel ordering (see camino/man/man1/camino.1). Here is the command: \\

Changed lines 85-105 from:

>> fid = fopen('PigBrainEIG.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> % Reshape into slices; 12 values per voxel \\ >> eigs = reshape(data, 12, 128,128,10); \\ >> % Elements 2, 3 and 4 are the principal eigenvector \\ >> fapds = abs(permute(eigs(2:4,:,:,:), [2,3,4,1])); \\ >> % Weight with fa. \\ >> fapds(:,:,:,1) = fapds(:,:,:,1).*fa; \\ >> fapds(:,:,:,2) = fapds(:,:,:,2).*fa; \\ >> fapds(:,:,:,3) = fapds(:,:,:,3).*fa; \\ >> fapds=permute(fapds, [1,2,4,3]); \\ >> % Transpose to standard orientation. \\ >> for i=1:3; for j=1:10 \\ fapds(:,:,i,j)=flipud(fapds(:,:,i,j)'); \\ end; end \\ >> imshow(fapds(:,:,:,5)) \\ >> for i=1:10 \\ imwrite(fapds(:,:,:,i), sprintf('pds%02i.png', i)); \\ end \\

to:

>> fid = fopen('PigBrainEIG.Bdouble', 'r', 'b');
>> data = fread(fid, 'double'); \\ >> fclose(fid);
>> % Reshape into slices; 12 values per voxel
>> eigs = reshape(data, 12, 128,128,10);
>> % Elements 2, 3 and 4 are the principal eigenvector
>> fapds = abs(permute(eigs(2:4,:,:,:), [2,3,4,1]));
>> % Weight with fa.
>> fapds(:,:,:,1) = fapds(:,:,:,1).*fa;
>> fapds(:,:,:,2) = fapds(:,:,:,2).*fa;
>> fapds(:,:,:,3) = fapds(:,:,:,3).*fa;
>> fapds=permute(fapds, [1,2,4,3]);
>> % Transpose to standard orientation.
>> for i=1:3;
>> for j=1:10
>> fapds(:,:,i,j)=flipud(fapds(:,:,i,j)');
>> end;
>> end
>> imshow(fapds(:,:,:,5))
>> for i=1:10
>> imwrite(fapds(:,:,:,i), sprintf('pds%02i.png', i));
>> end\\

Changed line 127 from:

files/pig/PigBrain05_PAS.png \\

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/PigBrain05_PAS.png \\

Changed line 142 from:

files/pig/VC_ThreshSelect.png \\

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/VC_ThreshSelect.png \\

Changed line 147 from:

\\ voxelclassify -inputfile PigBrain.Bfloat -bgthresh 500 -schemefile PigBrain.scheme -order 4 -ftest 1.0E-20 1.0E-7 1.0E-7 > PigBrainVC.Bint \\

to:

voxelclassify -inputfile PigBrain.Bfloat -bgthresh 500 -schemefile PigBrain.scheme -order 4 -ftest 1.0E-20 1.0E-7 1.0E-7 > PigBrainVC.Bint \\

Changed line 149 from:

The output of this command has type int and every voxel contains either -1 for background or 0, 2 or 4 indicating the classification order. A nice visualization of the results highlights the fibre crossings on the fractional anisotropy map. This is simple in matlab: \\

to:

The output of this command has type int and every voxel contains either -1 for background or 0, 2 or 4 indicating the classification order. A nice visualization of the results highlights the fibre crossings on the fractional anisotropy map. This is simple in Matlab: \\

Changed line 152 from:

files/pig/favc05.png

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/favc05.png

October 22, 2009, at 02:08 AM by shahrum -

Changed line 58 from:

>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));\\

to:

>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));\\

Changed line 68 from:

>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));\\

to:

>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));\\

Changed line 72 from:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/trd05.pngfiles/pig/fa05.png \\

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/fa05.png \\

October 22, 2009, at 02:07 AM by shahrum -

Changed lines 54-55 from:

>> % Show one slice \\ >> imshow(flipud(trd(:,:,5)'));\\

to:

>> % Show one slice
>> imshow(flipud(trd(:,:,5)'));\\

Changed line 58 from:

>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));\\

to:

>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));\\

October 22, 2009, at 02:06 AM by shahrum -

Changed lines 50-53 from:

>> data = fread(fid, 'double');\\ >> fclose(fid);\\ >> % Normalize and reshape into slices\\ >> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data));\\

to:

>> data = fread(fid, 'double');
>> fclose(fid);
>> % Normalize and reshape into slices
>> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data));\\

Changed line 57 from:

>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));\\

to:

>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));\\

Changed line 61 from:

>> fid = fopen('PigBrainFA.Bdouble', 'r', 'b');\\

to:

>> fid = fopen('PigBrainFA.Bdouble', 'r', 'b');\\

Changed line 67 from:

>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));\\

to:

>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));\\

October 22, 2009, at 02:05 AM by shahrum -

Changed lines 50-63 from:

>> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> % Normalize and reshape into slices \\ >> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> % Show one slice \\ >> imshow(flipud(trd(:,:,5)')); \\ >> % Write out all slices as .png images \\ >> for i=1:10 \\ >> -> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i)); \\ >> end \\ >>
>> % Repeat for FA
>> fid = fopen('PigBrainFA.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double');
>> fclose(fid); \\

to:

>> data = fread(fid, 'double');\\ >> fclose(fid);\\ >> % Normalize and reshape into slices\\ >> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data));\\ >> % Show one slice \\ >> imshow(flipud(trd(:,:,5)'));
>> % Write out all slices as .png images
>> for i=1:10
>> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i));\\ >> end

>> % Repeat for FA
>> fid = fopen('PigBrainFA.Bdouble', 'r', 'b');\\ >> data = fread(fid, 'double');
>> fclose(fid);\\

Changed lines 65-68 from:

>> imshow(flipud(fa(:,:,5)')); \\ >> for i=1:10 \\

>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));
>> end \\

to:

>> imshow(flipud(fa(:,:,5)'));
>> for i=1:10
>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));
>> end\\

October 22, 2009, at 02:02 AM by shahrum -

Changed line 44 from:

trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble

to:

trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble\\

October 22, 2009, at 02:01 AM by shahrum -

Changed lines 44-45 from:

@@trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble@@

to:

trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble

October 22, 2009, at 02:01 AM by shahrum -

Changed lines 44-45 from:

trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble \\ fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble

to:

@@trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble@@

Changed line 49 from:

>> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b');\\

to:

>> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b');\\

October 22, 2009, at 01:59 AM by shahrum -

Changed line 48 from:

>> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b'); \\

to:

>> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b');\\

October 22, 2009, at 01:58 AM by shahrum -

Changed lines 11-13 from:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the format. The simple perl script camino/PointSets/PointSetToScheme converts pigPoints.txt to a scheme file:

PointSets/PointSetToScheme pigPoints.txt 375330 3 0.0223 \\ \\

to:

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the format. The simple perl script camino/PointSets/PointSetToScheme converts pigPoints.txt to a scheme file:

PointSets/PointSetToScheme pigPoints.txt 375330 3 0.0223

Changed lines 27-33 from:

To fit the diffusion tensor in each voxel we can use dtfit or modelfit: \\

to:

To fit the diffusion tensor in each voxel we can use dtfit or modelfit:

dtfit PigBrain.Bfloat PigBrain.scheme > PigBrainDT.Bdouble

or equivalently:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 > PigBrainDT.Bdouble \\

Changed lines 35-36 from:

dtfit PigBrain.Bfloat PigBrain.scheme > PigBrainDT.Bdouble \\
or equivalently: \\

to:

Both commands give the same output. Note that the default output data type is big-endian doubles. The modelfit program has more options (see the modelfit man page) than dtfit, such as thresholding the b=0 measurements to remove the noisy background you can see in many of the images below, like this: \\

Changed line 37 from:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 > PigBrainDT.Bdouble \\

to:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 -bgthresh 500 > PigBrainDT.Bdouble \\

Changed lines 39-43 from:

Both commands give the same output. Note that the default output data type is big-endian doubles. The modelfit program has more options (see the modelfit man page) than dtfit, such as thresholding the b=0 measurements to remove the noisy background you can see in many of the images below, like this:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 -bgthresh 500 > PigBrainDT.Bdouble

to:

Changed lines 42-43 from:

From the diffusion tensor, we can compute simple statistics such as Trace(D), fractional anisotropy and principal directions. The commands trd and fa compute the first two:
\\ trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble \\ fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble \\

to:

From the diffusion tensor, we can compute simple statistics such as Trace(D), fractional anisotropy and principal directions. The commands trd and fa compute the first two:

trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble \\ fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble

A simple way to look at the results is to load the data files into matlab and display with imshow (use image instead if you do not have the image-processing toolbox) or write out slice by slice:

>> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> % Normalize and reshape into slices \\ >> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> % Show one slice \\ >> imshow(flipud(trd(:,:,5)')); \\ >> % Write out all slices as .png images \\ >> for i=1:10 \\ >> -> imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i)); \\ >> end \\ >>
>> % Repeat for FA
>> fid = fopen('PigBrainFA.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double');
>> fclose(fid); \\ >> fa = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data));
>> imshow(flipud(fa(:,:,5)')); \\ >> for i=1:10 \\

>> imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i));
>> end \\

Deleted lines 68-70:

A simple way to look at the results is to load the data files into matlab and display with imshow (use image instead if you do not have the image-processing toolbox) or write out slice by slice:
\\ >> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> % Normalize and reshape into slices \\ >> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> % Show one slice \\ >> imshow(flipud(trd(:,:,5)')); \\ >> % Write out all slices as .png images \\ >> for i=1:10 \\ imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i)); \\ end \\ >> \\ >> % Repeat for FA \\ >> fid = fopen('PigBrainFA.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> fa = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> imshow(flipud(fa(:,:,5)')); \\ >> for i=1:10 \\ imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i)); \\ end
\\

Changed line 70 from:

files/pig/trd05.pngfiles/pig/fa05.png \\

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/trd05.pngfiles/pig/fa05.png \\

Changed lines 74-75 from:

analyzeheader -datadims 128 128 10 -voxeldims 0.5 0.5 0.5 -datatype double > fa.hdr \\ \\ cp PigBrainFA.Bdouble fa.img \\

to:

analyzeheader -datadims 128 128 10 -voxeldims 0.5 0.5 0.5 -datatype double > fa.hdr \\ cp PigBrainFA.Bdouble fa.img \\

Changed line 81 from:

We can generate colour-coded direction maps easily from the output in matlab: \\

to:

We can generate colour-coded direction maps easily from the output in Matlab: \\

Changed line 86 from:

files/pig/pds05.png

to:

http://www.cs.ucl.ac.uk/research/medic/camino/tutorials/files/pig/pds05.png

September 04, 2009, at 04:32 PM by shahrum -

Added lines 1-106:

Tutorial: detecting crossing fibres

This tutorial uses Camino to generate fractional-anisotropy, colour-coded direction and voxel-classification maps. The voxel classification maps indicate which voxels exhibit isotropic diffusion, anisotropic Gaussian diffusion (adequately modelled by a single diffusion tensor) and non-Gaussian displacements (the diffusion tensor model fails), which are likely fibre crossings.

Data

The data set in this case study comes from a post-mortem porcine brain and was kindly provided by Tim Dyrby from The Danish Research Centre for Magnetic Resonance, Copenhagen University Hospital, Hvidovre, Denmark. The data set is intended for demonstration and testing of the Camino software package only. If you wish to use this data for other purposes, such as publications, other web sites or commercial enterprises, please contact Tim at (timd|at|drcmr.dk) before doing so. You can download the data set . The acquisition uses a spherical acquisition scheme with 61 unique gradient directions. The gradient directions come from the file camino/PointSets/Elec061.txt and minimize the electrostatic energy for 61 pairs of equal-and-opposite point electrical charges free to move on the surface of the sphere. The procedure outlined in [Jansons and Alexander, Inverse Problems, Vol 19, pp 1031-1046, 2003] computed them. The diffusion-weighting gradient strength G = 61 mT m^-1; pulse onset separation DELTA = 30ms; pulse width delta = 23ms; TE = 60ms; b = 3146 s mm^-2; NEX=2. The protocol also includes 3 acquisitions with no diffusion weighting. The image has 10 slices with in-plane resolution 128x128; voxel size 0.5x0.5x0.5 mm³.

Scheme File

Before we can process the data, we must create a scheme file that tells Camino about the acquisition sequence. The scheme file is a text file; the Camino man page (camino/man/man1/camino.1) specifies the format. The simple perl script camino/PointSets/PointSetToScheme converts pigPoints.txt to a scheme file:

PointSets/PointSetToScheme pigPoints.txt 375330 3 0.0223 \\
The first argument is the point set file; the second is the wavenumber q = gamma*delta*G = 2.6751987E8 * 0.023 * 0.061, where gamma is the gyromagnetic ratio of protons in water; the third argument is the number of b=0 images and the last is the diffusion time. Here we use the standard approximation of the diffusion time for finite pulse widths t = DELTA - delta/3 = 0.03 - 0.23/3. This is the resulting . A scheme file specifies the gradient directions along with the zero directions. The zero directions should always be positioned before the non-zero directions.

Data File

Each acquisition from the scanner is in an analyze format file NEX2_s_20060127_001_semsdw_15_image??, where ?? goes from 01 to 64; 01-03 are the b=0 images and 04-64 are the diffusion-weighted images with gradients in the same order as the scheme file. Each file stores image intensities as little-endian 2-byte short integers. We must generate a data file in camino format, which requires switching to big endian representation and combining all the analyze files into a single file with voxel ordering (see camino/man/man1/camino.1). Here is the command:

cat *.img | shredder 0 -2 0 | scanner2voxel -voxels 163840 -components 64 -inputdatatype short > PigBrain.Bfloat

The shredder command switches the byte ordering and scanner2voxel ensures that all measurements for each voxel are contiguous in the data file. The datafile PigBrain.Bfloat contains big-endian 4-byte floating-point data.

Fit the diffusion tensor

To fit the diffusion tensor in each voxel we can use dtfit or modelfit:

dtfit PigBrain.Bfloat PigBrain.scheme > PigBrainDT.Bdouble \\
or equivalently:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 > PigBrainDT.Bdouble

Both commands give the same output. Note that the default output data type is big-endian doubles. The modelfit program has more options (see the modelfit man page) than dtfit, such as thresholding the b=0 measurements to remove the noisy background you can see in many of the images below, like this:

modelfit -inputfile PigBrain.Bfloat -schemefile PigBrain.scheme -inversion 1 -bgthresh 500 > PigBrainDT.Bdouble

Simple tensor statistics

From the diffusion tensor, we can compute simple statistics such as Trace(D), fractional anisotropy and principal directions. The commands trd and fa compute the first two:
\\ trd < PigBrainDT.Bdouble > PigBrainTrD.Bdouble \\ fa < PigBrainDT.Bdouble > PigBrainFA.Bdouble

A simple way to look at the results is to load the data files into matlab and display with imshow (use image instead if you do not have the image-processing toolbox) or write out slice by slice:
\\ >> fid = fopen('PigBrainTrD.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> % Normalize and reshape into slices \\ >> trd = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> % Show one slice \\ >> imshow(flipud(trd(:,:,5)')); \\ >> % Write out all slices as .png images \\ >> for i=1:10 \\ imwrite(flipud(trd(:,:,i)'), sprintf('trd%02i.png', i)); \\ end \\ >> \\ >> % Repeat for FA \\ >> fid = fopen('PigBrainFA.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> fa = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> imshow(flipud(fa(:,:,5)')); \\ >> for i=1:10 \\ imwrite(flipud(fa(:,:,i)'), sprintf('fa%02i.png', i)); \\ end

That outputs the following images (slice 5 only):
files/pig/trd05.pngfiles/pig/fa05.png

Alternatively, you can convert the camino data files to analyze format, just by adding a header, and view them in, for example, MRICro:

analyzeheader -datadims 128 128 10 -voxeldims 0.5 0.5 0.5 -datatype double > fa.hdr \\ \\ cp PigBrainFA.Bdouble fa.img

dteig computes the eigensystem of each tensor:

dteig < PigBrainDT.Bdouble > PigBrainEIG.Bdouble

We can generate colour-coded direction maps easily from the output in matlab:

>> fid = fopen('PigBrainEIG.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> % Reshape into slices; 12 values per voxel \\ >> eigs = reshape(data, 12, 128,128,10); \\ >> % Elements 2, 3 and 4 are the principal eigenvector \\ >> fapds = abs(permute(eigs(2:4,:,:,:), [2,3,4,1])); \\ >> % Weight with fa. \\ >> fapds(:,:,:,1) = fapds(:,:,:,1).*fa; \\ >> fapds(:,:,:,2) = fapds(:,:,:,2).*fa; \\ >> fapds(:,:,:,3) = fapds(:,:,:,3).*fa; \\ >> fapds=permute(fapds, [1,2,4,3]); \\ >> % Transpose to standard orientation. \\ >> for i=1:3; for j=1:10 \\ fapds(:,:,i,j)=flipud(fapds(:,:,i,j)'); \\ end; end \\ >> imshow(fapds(:,:,:,5)) \\ >> for i=1:10 \\ imwrite(fapds(:,:,:,i), sprintf('pds%02i.png', i)); \\ end

That outputs the following image (slice 5 only):
files/pig/pds05.png

Diffusion tensor plot

The command sfplot creates plots of diffusion MRI reconstructions over slices. In particular, it can create plots of the diffusion tensor over a slice. It can also create plots of other reconstructed features, such as the PAS or q-ball ODF; see for example the Parallel PASMRI case study or the sfplot man page. The following commands create a plot of the DT over slice 5 of the pig brain data set. We start by creating a fractional anisotropy map to use as background for the plot:

fa < PigBrainDT.Bdouble | shredder $((128*128*8*4)) $((128*128*8)) $((128*128*8*10)) > PigBrainFA_SLICE05.Bdouble

The command above uses shredder to pick out the fifth slice of the FA map over the whole volume; each slice contains 128*128*8 bytes, as the FA in each of the 128*128 voxels per slice is an 8-byte double. Now create the plot itself:

shredder $((128*128*8*8*4)) $((128*128*8*8)) $((128*128*8*8*10)) < PigBrainDT.Bdouble | sfplot -xsize 128 -ysize 128 -inputmodel dt -backdrop PigBrainFA_SLICE05.Bdouble -dircolcode -maxd 0.4E-9 -projection 2 1 > PigBrain05_PAS.rgb

The man page for sfplot explains the options in detail. Briefly, the first two (-xsize 128 -ysize 128) define the size of the image array; -inputmodel dt says that the program is reading in diffusion tensor data; -backdrop PigBrainFA_SLICE05.Bdouble says to use the FA map we just created as the background; -dircolcode defines the colour scheme to use and encodes directions as colours so that left-right oriented tensors are predominantly red, anterior-posterior green and superior inferior blue (other schemes are available); -maxd 0.4E-9 specifies the scale of the tensors, as the default is larger as it was chosen for in-vivo imaging; -projection 2 1 specifies the orientations of the DTs with respect to the image. In practice, you need to experiment with some of these settings to get good images for particular acquisitions. The shredder command above picks out the fifth slice of the diffusion tensor data set.

The command outputs a figure in raw rgb format. You can convert it to something more convenient using, for example, ImageMagick's convert command:

convert -size 2048x2048 -depth 8 PigBrain05_PAS.rgb PigBrain05_PAS.png

Here is the figure output (click to see the full size image):
files/pig/PigBrain05_PAS.png

To make the figure look a bit nicer, I've used a version of the PigBrainDT.Bdouble file with the background threshold set to 500.

Voxel Classification

The command voxelclassify classifies each voxel as isotropic, anisotropic Gaussian or non-Gaussian using the algorithm in [Alexander, Barker and Arridge, MRM, Vol 48, pp 331-340, 2002]. In isotropic and anisotropic Gaussian voxels, the diffusion tensor model works well. Non-Gaussian voxels are likely fibre crossings where the diffusion tensor fails. The algorithm requires you to select thresholds for separating the three classes. First run:

voxelclassify -inputfile PigBrain.Bfloat -bgthresh 500 -schemefile PigBrain.scheme -order 4 > PigBrainVC.Bdouble

To choose thresholds, run the vcthreshselect program giving it the output of the command above:

vcthreshselect -inputfile PigBrainVC.Bdouble -datadims 128 128 10 -order 4

A window appears that looks like this:
files/pig/VC_ThreshSelect.png

Black voxels are isotropic ("Order 0") or background, dark gray are anisotropic Gaussian ("Order 2") and light gray are non-Gaussian ("Order 4"). The light gray ones are likely fibre crossings. Scroll through the slices using the list on the right. Set the thresholds using the scroll bars or text fields at the bottom. At good threshold settings, the order 4 voxels should be fairly clustered, ie as few isolated order 4 voxels as possible with still a reasonable fraction of order 4 voxels overall. Around 5 to 10 percent order 4 voxels is typically about right.

Once you have chosen thresholds, make a note of them and generate the classified volume using voxelclassify with specified thresholds:
\\ voxelclassify -inputfile PigBrain.Bfloat -bgthresh 500 -schemefile PigBrain.scheme -order 4 -ftest 1.0E-20 1.0E-7 1.0E-7 > PigBrainVC.Bint

The output of this command has type int and every voxel contains either -1 for background or 0, 2 or 4 indicating the classification order. A nice visualization of the results highlights the fibre crossings on the fractional anisotropy map. This is simple in matlab:

>> fid = fopen('PigBrainVC.Bint', 'r', 'b'); \\ >> vc=fread(fid, 'int'); \\ >> fclose(fid); \\ >> vc=reshape(vc, 128,128,10); \\ >> % Repeat for FA \\ >> fid = fopen('PigBrainFA.Bdouble', 'r', 'b'); \\ >> data = fread(fid, 'double'); \\ >> fclose(fid); \\ >> fa = (reshape(data, 128,128,10) - min(data))/(max(data) - min(data)); \\ >> for i=1:10 \\ im(:,:,1) = fa(:,:,i).*(1-0.5*(vc(:,:,i)>=4)) + 0.5*(vc(:,:,i)>=4); \\ im(:,:,2) = fa(:,:,i).*(1-0.5*(vc(:,:,i)>=4)) + 0.5*(vc(:,:,i)>=4); \\ im(:,:,3) = fa(:,:,i); \\ imwrite(im, sprintf('/tmp/favc%02i.png', i)); \\ end \\\\ Here is the image for slice 5:
files/pig/favc05.png

Tutorial: detecting crossing fibres

Tutorials.PigBrainTutorial History

Diffusion tensor plot

Data File

Simple tensor statistics

Tutorial: detecting crossing fibres

Tutorial: detecting crossing fibres

Data

Scheme File

Data File

Fit the diffusion tensor

Simple tensor statistics

Diffusion tensor plot

Voxel Classification