Skip to content

GitLab

Commit ce12fa61 authored by Iyad Bagari's avatar Iyad Bagari
Browse files

Initial commit

parents
Hide whitespace changes
Inline Side-by-side
Showing with 3194 additions and 0 deletions
.DS_Store 0 → 100644
View file @ ce12fa61
File added
NIfTI_20140122/FAQ.pdf 0 → 100644
View file @ ce12fa61
File added
NIfTI_20140122/NIfTI_tools.pdf 0 → 100644
View file @ ce12fa61
File added
NIfTI_20140122/UseANALYZE.pdf 0 → 100644
View file @ ce12fa61
File added
NIfTI_20140122/affine.m 0 → 100644
View file @ ce12fa61
% Using 2D or 3D affine matrix to rotate, translate, scale, reflect and
% shear a 2D image or 3D volume. 2D image is represented by a 2D matrix,
% 3D volume is represented by a 3D matrix, and data type can be real
% integer or floating-point.
%
% You may notice that MATLAB has a function called 'imtransform.m' for
% 2D spatial transformation. However, keep in mind that 'imtransform.m'
% assumes y for the 1st dimension, and x for the 2nd dimension. They are
% equivalent otherwise.
%
% In addition, if you adjust the 'new_elem_size' parameter, this 'affine.m'
% is equivalent to 'interp2.m' for 2D image, and equivalent to 'interp3.m'
% for 3D volume.
%
% Usage: [new_img new_M] = ...
% affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);
%
% old_img - original 2D image or 3D volume. We assume x for the 1st
% dimension, y for the 2nd dimension, and z for the 3rd
% dimension.
%
% old_M - a 3x3 2D affine matrix for 2D image, or a 4x4 3D affine
% matrix for 3D volume. We assume x for the 1st dimension,
% y for the 2nd dimension, and z for the 3rd dimension.
%
% new_elem_size (optional) - size of voxel along x y z direction for
% a transformed 3D volume, or size of pixel along x y for
% a transformed 2D image. We assume x for the 1st dimension
% y for the 2nd dimension, and z for the 3rd dimension.
% 'new_elem_size' is 1 if it is default or empty.
%
% You can increase its value to decrease the resampling rate,
% and make the 2D image or 3D volume more coarse. It works
% just like 'interp3'.
%
% verbose (optional) - 1, 0
% 1: show transforming progress in percentage
% 2: progress will not be displayed
% 'verbose' is 1 if it is default or empty.
%
% bg (optional) - background voxel intensity in any extra corner that
% is caused by the interpolation. 0 in most cases. If it is
% default or empty, 'bg' will be the average of two corner
% voxel intensities in original data.
%
% method (optional) - 1, 2, or 3
% 1: for Trilinear interpolation
% 2: for Nearest Neighbor interpolation
% 3: for Fischer's Bresenham interpolation
% 'method' is 1 if it is default or empty.
%
% new_img - transformed 2D image or 3D volume
%
% new_M - transformed affine matrix
%
% Example 1 (3D rotation):
% load mri.mat; old_img = double(squeeze(D));
% old_M = [0.88 0.5 3 -90; -0.5 0.88 3 -126; 0 0 2 -72; 0 0 0 1];
% new_img = affine(old_img, old_M, 2);
% [x y z] = meshgrid(1:128,1:128,1:27);
% sz = size(new_img);
% [x1 y1 z1] = meshgrid(1:sz(2),1:sz(1),1:sz(3));
% figure; slice(x, y, z, old_img, 64, 64, 13.5);
% shading flat; colormap(map); view(-66, 66);
% figure; slice(x1, y1, z1, new_img, sz(1)/2, sz(2)/2, sz(3)/2);
% shading flat; colormap(map); view(-66, 66);
%
% Example 2 (2D interpolation):
% load mri.mat; old_img=D(:,:,1,13)';
% old_M = [1 0 0; 0 1 0; 0 0 1];
% new_img = affine(old_img, old_M, [.2 .4]);
% figure; image(old_img); colormap(map);
% figure; image(new_img); colormap(map);
%
% This program is inspired by:
% SPM5 Software from Wellcome Trust Centre for Neuroimaging
% http://www.fil.ion.ucl.ac.uk/spm/software
% Fischer, J., A. del Rio (2004). A Fast Method for Applying Rigid
% Transformations to Volume Data, WSCG2004 Conference.
% http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function [new_img, new_M] = affine(old_img, old_M, new_elem_size, verbose, bg, method)
if ~exist('old_img','var') | ~exist('old_M','var')
error('Usage: [new_img new_M] = affine(old_img, old_M, [new_elem_size], [verbose], [bg], [method]);');
end
if ndims(old_img) == 3
if ~isequal(size(old_M),[4 4])
error('old_M should be a 4x4 affine matrix for 3D volume.');
end
elseif ndims(old_img) == 2
if ~isequal(size(old_M),[3 3])
error('old_M should be a 3x3 affine matrix for 2D image.');
end
else
error('old_img should be either 2D image or 3D volume.');
end
if ~exist('new_elem_size','var') | isempty(new_elem_size)
new_elem_size = [1 1 1];
elseif length(new_elem_size) < 2
new_elem_size = new_elem_size(1)*ones(1,3);
elseif length(new_elem_size) < 3
new_elem_size = [new_elem_size(:); 1]';
end
if ~exist('method','var') | isempty(method)
method = 1;
elseif ~exist('bresenham_line3d.m','file') & method == 3
error([char(10) char(10) 'Please download 3D Bresenham''s line generation program from:' char(10) char(10) 'http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21057' char(10) char(10) 'to test Fischer''s Bresenham interpolation method.' char(10) char(10)]);
end
% Make compatible to MATLAB earlier than version 7 (R14), which
% can only perform arithmetic on double data type
%
old_img = double(old_img);
old_dim = size(old_img);
if ~exist('bg','var') | isempty(bg)
bg = mean([old_img(1) old_img(end)]);
end
if ~exist('verbose','var') | isempty(verbose)
verbose = 1;
end
if ndims(old_img) == 2
old_dim(3) = 1;
old_M = old_M(:, [1 2 3 3]);
old_M = old_M([1 2 3 3], :);
old_M(3,:) = [0 0 1 0];
old_M(:,3) = [0 0 1 0]';
end
% Vertices of img in voxel
%
XYZvox = [ 1 1 1
1 1 old_dim(3)
1 old_dim(2) 1
1 old_dim(2) old_dim(3)
old_dim(1) 1 1
old_dim(1) 1 old_dim(3)
old_dim(1) old_dim(2) 1
old_dim(1) old_dim(2) old_dim(3) ]';
old_R = old_M(1:3,1:3);
old_T = old_M(1:3,4);
% Vertices of img in millimeter
%
XYZmm = old_R*(XYZvox-1) + repmat(old_T, [1, 8]);
% Make scale of new_M according to new_elem_size
%
new_M = diag([new_elem_size 1]);
% Make translation so minimum vertex is moved to [1,1,1]
%
new_M(1:3,4) = round( min(XYZmm,[],2) );
% New dimensions will be the maximum vertices in XYZ direction (dim_vox)
% i.e. compute dim_vox via dim_mm = R*(dim_vox-1)+T
% where, dim_mm = round(max(XYZmm,[],2));
%
new_dim = ceil(new_M(1:3,1:3) \ ( round(max(XYZmm,[],2))-new_M(1:3,4) )+1)';
% Initialize new_img with new_dim
%
new_img = zeros(new_dim(1:3));
% Mask out any changes from Z axis of transformed volume, since we
% will traverse it voxel by voxel below. We will only apply unit
% increment of mask_Z(3,4) to simulate the cursor movement
%
% i.e. we will use mask_Z * new_XYZvox to replace new_XYZvox
%
mask_Z = diag(ones(1,4));
mask_Z(3,3) = 0;
% It will be easier to do the interpolation if we invert the process
% by not traversing the original volume. Instead, we traverse the
% transformed volume, and backproject each voxel in the transformed
% volume back into the original volume. If the backprojected voxel
% in original volume is within its boundary, the intensity of that
% voxel can be used by the cursor location in the transformed volume.
%
% First, we traverse along Z axis of transformed volume voxel by voxel
%
for z = 1:new_dim(3)
if verbose & ~mod(z,10)
fprintf('%.2f percent is done.\n', 100*z/new_dim(3));
end
% We need to find out the mapping from voxel in the transformed
% volume (new_XYZvox) to voxel in the original volume (old_XYZvox)
%
% The following equation works, because they all equal to XYZmm:
% new_R*(new_XYZvox-1) + new_T == old_R*(old_XYZvox-1) + old_T
%
% We can use modified new_M1 & old_M1 to substitute new_M & old_M
% new_M1 * new_XYZvox == old_M1 * old_XYZvox
%
% where: M1 = M; M1(:,4) = M(:,4) - sum(M(:,1:3),2);
% and: M(:,4) == [T; 1] == sum(M1,2)
%
% Therefore: old_XYZvox = old_M1 \ new_M1 * new_XYZvox;
%
% Since we are traverse Z axis, and new_XYZvox is replaced
% by mask_Z * new_XYZvox, the above formula can be rewritten
% as: old_XYZvox = old_M1 \ new_M1 * mask_Z * new_XYZvox;
%
% i.e. we find the mapping from new_XYZvox to old_XYZvox:
% M = old_M1 \ new_M1 * mask_Z;
%
% First, compute modified old_M1 & new_M1
%
old_M1 = old_M; old_M1(:,4) = old_M(:,4) - sum(old_M(:,1:3),2);
new_M1 = new_M; new_M1(:,4) = new_M(:,4) - sum(new_M(:,1:3),2);
% Then, apply unit increment of mask_Z(3,4) to simulate the
% cursor movement
%
mask_Z(3,4) = z;
% Here is the mapping from new_XYZvox to old_XYZvox
%
M = old_M1 \ new_M1 * mask_Z;
switch method
case 1
new_img(:,:,z) = trilinear(old_img, new_dim, old_dim, M, bg);
case 2
new_img(:,:,z) = nearest_neighbor(old_img, new_dim, old_dim, M, bg);
case 3
new_img(:,:,z) = bresenham(old_img, new_dim, old_dim, M, bg);
end
end; % for z
if ndims(old_img) == 2
new_M(3,:) = [];
new_M(:,3) = [];
end
return; % affine
%--------------------------------------------------------------------
function img_slice = trilinear(img, dim1, dim2, M, bg)
img_slice = zeros(dim1(1:2));
TINY = 5e-2; % tolerance
% Dimension of transformed 3D volume
%
xdim1 = dim1(1);
ydim1 = dim1(2);
% Dimension of original 3D volume
%
xdim2 = dim2(1);
ydim2 = dim2(2);
zdim2 = dim2(3);
% initialize new_Y accumulation
%
Y2X = 0;
Y2Y = 0;
Y2Z = 0;
for y = 1:ydim1
% increment of new_Y accumulation
%
Y2X = Y2X + M(1,2); % new_Y to old_X
Y2Y = Y2Y + M(2,2); % new_Y to old_Y
Y2Z = Y2Z + M(3,2); % new_Y to old_Z
% backproject new_Y accumulation and translation to old_XYZ
%
old_X = Y2X + M(1,4);
old_Y = Y2Y + M(2,4);
old_Z = Y2Z + M(3,4);
for x = 1:xdim1
% accumulate the increment of new_X, and apply it
% to the backprojected old_XYZ
%
old_X = M(1,1) + old_X ;
old_Y = M(2,1) + old_Y ;
old_Z = M(3,1) + old_Z ;
% within boundary of original image
%
if ( old_X > 1-TINY & old_X < xdim2+TINY & ...
old_Y > 1-TINY & old_Y < ydim2+TINY & ...
old_Z > 1-TINY & old_Z < zdim2+TINY )
% Calculate distance of old_XYZ to its neighbors for
% weighted intensity average
%
dx = old_X - floor(old_X);
dy = old_Y - floor(old_Y);
dz = old_Z - floor(old_Z);
x000 = floor(old_X);
x100 = x000 + 1;
if floor(old_X) < 1
x000 = 1;
x100 = x000;
elseif floor(old_X) > xdim2-1
x000 = xdim2;
x100 = x000;
end
x010 = x000;
x001 = x000;
x011 = x000;
x110 = x100;
x101 = x100;
x111 = x100;
y000 = floor(old_Y);
y010 = y000 + 1;
if floor(old_Y) < 1
y000 = 1;
y100 = y000;
elseif floor(old_Y) > ydim2-1
y000 = ydim2;
y010 = y000;
end
y100 = y000;
y001 = y000;
y101 = y000;
y110 = y010;
y011 = y010;
y111 = y010;
z000 = floor(old_Z);
z001 = z000 + 1;
if floor(old_Z) < 1
z000 = 1;
z001 = z000;
elseif floor(old_Z) > zdim2-1
z000 = zdim2;
z001 = z000;
end
z100 = z000;
z010 = z000;
z110 = z000;
z101 = z001;
z011 = z001;
z111 = z001;
x010 = x000;
x001 = x000;
x011 = x000;
x110 = x100;
x101 = x100;
x111 = x100;
v000 = double(img(x000, y000, z000));
v010 = double(img(x010, y010, z010));
v001 = double(img(x001, y001, z001));
v011 = double(img(x011, y011, z011));
v100 = double(img(x100, y100, z100));
v110 = double(img(x110, y110, z110));
v101 = double(img(x101, y101, z101));
v111 = double(img(x111, y111, z111));
img_slice(x,y) = v000*(1-dx)*(1-dy)*(1-dz) + ...
v010*(1-dx)*dy*(1-dz) + ...
v001*(1-dx)*(1-dy)*dz + ...
v011*(1-dx)*dy*dz + ...
v100*dx*(1-dy)*(1-dz) + ...
v110*dx*dy*(1-dz) + ...
v101*dx*(1-dy)*dz + ...
v111*dx*dy*dz;
else
img_slice(x,y) = bg;
end % if boundary
end % for x
end % for y
return; % trilinear
%--------------------------------------------------------------------
function img_slice = nearest_neighbor(img, dim1, dim2, M, bg)
img_slice = zeros(dim1(1:2));
% Dimension of transformed 3D volume
%
xdim1 = dim1(1);
ydim1 = dim1(2);
% Dimension of original 3D volume
%
xdim2 = dim2(1);
ydim2 = dim2(2);
zdim2 = dim2(3);
% initialize new_Y accumulation
%
Y2X = 0;
Y2Y = 0;
Y2Z = 0;
for y = 1:ydim1
% increment of new_Y accumulation
%
Y2X = Y2X + M(1,2); % new_Y to old_X
Y2Y = Y2Y + M(2,2); % new_Y to old_Y
Y2Z = Y2Z + M(3,2); % new_Y to old_Z
% backproject new_Y accumulation and translation to old_XYZ
%
old_X = Y2X + M(1,4);
old_Y = Y2Y + M(2,4);
old_Z = Y2Z + M(3,4);
for x = 1:xdim1
% accumulate the increment of new_X and apply it
% to the backprojected old_XYZ
%
old_X = M(1,1) + old_X ;
old_Y = M(2,1) + old_Y ;
old_Z = M(3,1) + old_Z ;
xi = round(old_X);
yi = round(old_Y);
zi = round(old_Z);
% within boundary of original image
%
if ( xi >= 1 & xi <= xdim2 & ...
yi >= 1 & yi <= ydim2 & ...
zi >= 1 & zi <= zdim2 )
img_slice(x,y) = img(xi,yi,zi);
else
img_slice(x,y) = bg;
end % if boundary
end % for x
end % for y
return; % nearest_neighbor
%--------------------------------------------------------------------
function img_slice = bresenham(img, dim1, dim2, M, bg)
img_slice = zeros(dim1(1:2));
% Dimension of transformed 3D volume
%
xdim1 = dim1(1);
ydim1 = dim1(2);
% Dimension of original 3D volume
%
xdim2 = dim2(1);
ydim2 = dim2(2);
zdim2 = dim2(3);
for y = 1:ydim1
start_old_XYZ = round(M*[0 y 0 1]');
end_old_XYZ = round(M*[xdim1 y 0 1]');
[X Y Z] = bresenham_line3d(start_old_XYZ, end_old_XYZ);
% line error correction
%
% del = end_old_XYZ - start_old_XYZ;
% del_dom = max(del);
% idx_dom = find(del==del_dom);
% idx_dom = idx_dom(1);
% idx_other = [1 2 3];
% idx_other(idx_dom) = [];
%del_x1 = del(idx_other(1));
% del_x2 = del(idx_other(2));
% line_slope = sqrt((del_x1/del_dom)^2 + (del_x2/del_dom)^2 + 1);
% line_error = line_slope - 1;
% line error correction removed because it is too slow
for x = 1:xdim1
% rescale ratio
%
i = round(x * length(X) / xdim1);
if i < 1
i = 1;
elseif i > length(X)
i = length(X);
end
xi = X(i);
yi = Y(i);
zi = Z(i);
% within boundary of the old XYZ space
%
if ( xi >= 1 & xi <= xdim2 & ...
yi >= 1 & yi <= ydim2 & ...
zi >= 1 & zi <= zdim2 )
img_slice(x,y) = img(xi,yi,zi);
% if line_error > 1
% x = x + 1;
% if x <= xdim1
% img_slice(x,y) = img(xi,yi,zi);
% line_error = line_slope - 1;
% end
% end % if line_error
% line error correction removed because it is too slow
else
img_slice(x,y) = bg;
end % if boundary
end % for x
end % for y
return; % bresenham
NIfTI_20140122/bipolar.m 0 → 100644
View file @ ce12fa61
%BIPOLAR returns an M-by-3 matrix containing a blue-red colormap, in
% in which red stands for positive, blue stands for negative,
% and white stands for 0.
%
% Usage: cmap = bipolar(M, lo, hi, contrast); or cmap = bipolar;
%
% cmap: output M-by-3 matrix for BIPOLAR colormap.
% M: number of shades in the colormap. By default, it is the
% same length as the current colormap.
% lo: the lowest value to represent.
% hi: the highest value to represent.
%
% Inspired from the LORETA PASCAL program:
% http://www.unizh.ch/keyinst/NewLORETA
%
% jimmy@rotman-baycrest.on.ca
%
%----------------------------------------------------------------
function cmap = bipolar(M, lo, hi, contrast)
if ~exist('contrast','var')
contrast = 128;
end
if ~exist('lo','var')
lo = -1;
end
if ~exist('hi','var')
hi = 1;
end
if ~exist('M','var')
cmap = colormap;
M = size(cmap,1);
end
steepness = 10 ^ (1 - (contrast-1)/127);
pos_infs = 1e-99;
neg_infs = -1e-99;
doubleredc = [];
doublebluec = [];
if lo >= 0 % all positive
if lo == 0
lo = pos_infs;
end
for i=linspace(hi/M, hi, M)
t = exp(log(i/hi)*steepness);
doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
end
cmap = doubleredc;
elseif hi <= 0 % all negative
if hi == 0
hi = neg_infs;
end
for i=linspace(abs(lo)/M, abs(lo), M)
t = exp(log(i/abs(lo))*steepness);
doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
end
cmap = flipud(doublebluec);
else
if hi > abs(lo)
maxc = hi;
else
maxc = abs(lo);
end
for i=linspace(maxc/M, hi, round(M*hi/(hi-lo)))
t = exp(log(i/maxc)*steepness);
doubleredc = [doubleredc; [(1-t)+t,(1-t)+0,(1-t)+0]];
end
for i=linspace(maxc/M, abs(lo), round(M*abs(lo)/(hi-lo)))
t = exp(log(i/maxc)*steepness);
doublebluec = [doublebluec; [(1-t)+0,(1-t)+0,(1-t)+t]];
end
cmap = [flipud(doublebluec); doubleredc];
end
return; % bipolar
NIfTI_20140122/bresenham_line3d.m 0 → 100644
View file @ ce12fa61
% Generate X Y Z coordinates of a 3D Bresenham's line between
% two given points.
%
% A very useful application of this algorithm can be found in the
% implementation of Fischer's Bresenham interpolation method in my
% another program that can rotate three dimensional image volume
% with an affine matrix:
% http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=21080
%
% Usage: [X Y Z] = bresenham_line3d(P1, P2, [precision]);
%
% P1 - vector for Point1, where P1 = [x1 y1 z1]
%
% P2 - vector for Point2, where P2 = [x2 y2 z2]
%
% precision (optional) - Although according to Bresenham's line
% algorithm, point coordinates x1 y1 z1 and x2 y2 z2 should
% be integer numbers, this program extends its limit to all
% real numbers. If any of them are floating numbers, you
% should specify how many digits of decimal that you would
% like to preserve. Be aware that the length of output X Y
% Z coordinates will increase in 10 times for each decimal
% digit that you want to preserve. By default, the precision
% is 0, which means that they will be rounded to the nearest
% integer.
%
% X - a set of x coordinates on Bresenham's line
%
% Y - a set of y coordinates on Bresenham's line
%
% Z - a set of z coordinates on Bresenham's line
%
% Therefore, all points in XYZ set (i.e. P(i) = [X(i) Y(i) Z(i)])
% will constitute the Bresenham's line between P1 and P1.
%
% Example:
% P1 = [12 37 6]; P2 = [46 3 35];
% [X Y Z] = bresenham_line3d(P1, P2);
% figure; plot3(X,Y,Z,'s','markerface','b');
%
% This program is ported to MATLAB from:
%
% B.Pendleton. line3d - 3D Bresenham's (a 3D line drawing algorithm)
% ftp://ftp.isc.org/pub/usenet/comp.sources.unix/volume26/line3d, 1992
%
% Which is also referenced by:
%
% Fischer, J., A. del Rio (2004). A Fast Method for Applying Rigid
% Transformations to Volume Data, WSCG2004 Conference.
% http://wscg.zcu.cz/wscg2004/Papers_2004_Short/M19.pdf
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function [X,Y,Z] = bresenham_line3d(P1, P2, precision)
if ~exist('precision','var') | isempty(precision) | round(precision) == 0
precision = 0;
P1 = round(P1);
P2 = round(P2);
else
precision = round(precision);
P1 = round(P1*(10^precision));
P2 = round(P2*(10^precision));
end
d = max(abs(P2-P1)+1);
X = zeros(1, d);
Y = zeros(1, d);
Z = zeros(1, d);
x1 = P1(1);
y1 = P1(2);
z1 = P1(3);
x2 = P2(1);
y2 = P2(2);
z2 = P2(3);
dx = x2 - x1;
dy = y2 - y1;
dz = z2 - z1;
ax = abs(dx)*2;
ay = abs(dy)*2;
az = abs(dz)*2;
sx = sign(dx);
sy = sign(dy);
sz = sign(dz);
x = x1;
y = y1;
z = z1;
idx = 1;
if(ax>=max(ay,az)) % x dominant
yd = ay - ax/2;
zd = az - ax/2;
while(1)
X(idx) = x;
Y(idx) = y;
Z(idx) = z;
idx = idx + 1;
if(x == x2) % end
break;
end
if(yd >= 0) % move along y
y = y + sy;
yd = yd - ax;
end
if(zd >= 0) % move along z
z = z + sz;
zd = zd - ax;
end
x = x + sx; % move along x
yd = yd + ay;
zd = zd + az;
end
elseif(ay>=max(ax,az)) % y dominant
xd = ax - ay/2;
zd = az - ay/2;
while(1)
X(idx) = x;
Y(idx) = y;
Z(idx) = z;
idx = idx + 1;
if(y == y2) % end
break;
end
if(xd >= 0) % move along x
x = x + sx;
xd = xd - ay;
end
if(zd >= 0) % move along z
z = z + sz;
zd = zd - ay;
end
y = y + sy; % move along y
xd = xd + ax;
zd = zd + az;
end
elseif(az>=max(ax,ay)) % z dominant
xd = ax - az/2;
yd = ay - az/2;
while(1)
X(idx) = x;
Y(idx) = y;
Z(idx) = z;
idx = idx + 1;
if(z == z2) % end
break;
end
if(xd >= 0) % move along x
x = x + sx;
xd = xd - az;
end
if(yd >= 0) % move along y
y = y + sy;
yd = yd - az;
end
z = z + sz; % move along z
xd = xd + ax;
yd = yd + ay;
end
end
if precision ~= 0
X = X/(10^precision);
Y = Y/(10^precision);
Z = Z/(10^precision);
end
return; % bresenham_line3d
NIfTI_20140122/clip_nii.m 0 → 100644
View file @ ce12fa61
% CLIP_NII: Clip the NIfTI volume from any of the 6 sides
%
% Usage: nii = clip_nii(nii, [option])
%
% Inputs:
%
% nii - NIfTI volume.
%
% option - struct instructing how many voxel to be cut from which side.
%
% option.cut_from_L = ( number of voxel )
% option.cut_from_R = ( number of voxel )
% option.cut_from_P = ( number of voxel )
% option.cut_from_A = ( number of voxel )
% option.cut_from_I = ( number of voxel )
% option.cut_from_S = ( number of voxel )
%
% Options description in detail:
% ==============================
%
% cut_from_L: Number of voxels from Left side will be clipped.
%
% cut_from_R: Number of voxels from Right side will be clipped.
%
% cut_from_P: Number of voxels from Posterior side will be clipped.
%
% cut_from_A: Number of voxels from Anterior side will be clipped.
%
% cut_from_I: Number of voxels from Inferior side will be clipped.
%
% cut_from_S: Number of voxels from Superior side will be clipped.
%
% NIfTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function nii = clip_nii(nii, opt)
dims = abs(nii.hdr.dime.dim(2:4));
origin = abs(nii.hdr.hist.originator(1:3));
if isempty(origin) | all(origin == 0) % according to SPM
origin = round((dims+1)/2);
end
cut_from_L = 0;
cut_from_R = 0;
cut_from_P = 0;
cut_from_A = 0;
cut_from_I = 0;
cut_from_S = 0;
if nargin > 1 & ~isempty(opt)
if ~isstruct(opt)
error('option argument should be a struct');
end
if isfield(opt,'cut_from_L')
cut_from_L = round(opt.cut_from_L);
if cut_from_L >= origin(1) | cut_from_L < 0
error('cut_from_L cannot be negative or cut beyond originator');
end
end
if isfield(opt,'cut_from_P')
cut_from_P = round(opt.cut_from_P);
if cut_from_P >= origin(2) | cut_from_P < 0
error('cut_from_P cannot be negative or cut beyond originator');
end
end
if isfield(opt,'cut_from_I')
cut_from_I = round(opt.cut_from_I);
if cut_from_I >= origin(3) | cut_from_I < 0
error('cut_from_I cannot be negative or cut beyond originator');
end
end
if isfield(opt,'cut_from_R')
cut_from_R = round(opt.cut_from_R);
if cut_from_R > dims(1)-origin(1) | cut_from_R < 0
error('cut_from_R cannot be negative or cut beyond originator');
end
end
if isfield(opt,'cut_from_A')
cut_from_A = round(opt.cut_from_A);
if cut_from_A > dims(2)-origin(2) | cut_from_A < 0
error('cut_from_A cannot be negative or cut beyond originator');
end
end
if isfield(opt,'cut_from_S')
cut_from_S = round(opt.cut_from_S);
if cut_from_S > dims(3)-origin(3) | cut_from_S < 0
error('cut_from_S cannot be negative or cut beyond originator');
end
end
end
nii = make_nii(nii.img( (cut_from_L+1) : (dims(1)-cut_from_R), ...
(cut_from_P+1) : (dims(2)-cut_from_A), ...
(cut_from_I+1) : (dims(3)-cut_from_S), ...
:,:,:,:,:), nii.hdr.dime.pixdim(2:4), ...
[origin(1)-cut_from_L origin(2)-cut_from_P origin(3)-cut_from_I], ...
nii.hdr.dime.datatype, nii.hdr.hist.descrip);
return;
NIfTI_20140122/collapse_nii_scan.m 0 → 100644
View file @ ce12fa61
% Collapse multiple single-scan NIFTI files into a multiple-scan NIFTI file
%
% Usage: collapse_nii_scan(scan_file_pattern, [collapsed_fileprefix], [scan_file_folder])
%
% Here, scan_file_pattern should look like: 'myscan_0*.img'
% If collapsed_fileprefix is omit, 'multi_scan' will be used
% If scan_file_folder is omit, current file folder will be used
%
% The order of volumes in the collapsed file will be the order of
% corresponding filenames for those selected scan files.
%
% NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function collapse_nii_scan(scan_pattern, fileprefix, scan_path)
if ~exist('fileprefix','var')
fileprefix = 'multi_scan';
else
[tmp fileprefix] = fileparts(fileprefix);
end
if ~exist('scan_path','var'), scan_path = pwd; end
pnfn = fullfile(scan_path, scan_pattern);
file_lst = dir(pnfn);
flist = {file_lst.name};
flist = flist(:);
filename = flist{1};
v = version;
% Check file extension. If .gz, unpack it into temp folder
%
if length(filename) > 2 & strcmp(filename(end-2:end), '.gz')
if ~strcmp(filename(end-6:end), '.img.gz') & ...
~strcmp(filename(end-6:end), '.hdr.gz') & ...
~strcmp(filename(end-6:end), '.nii.gz')
error('Please check filename.');
end
if str2num(v(1:3)) < 7.1 | ~usejava('jvm')
error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.');
else
gzFile = 1;
end
else
if ~strcmp(filename(end-3:end), '.img') & ...
~strcmp(filename(end-3:end), '.hdr') & ...
~strcmp(filename(end-3:end), '.nii')
error('Please check filename.');
end
end
nii = load_untouch_nii(fullfile(scan_path,filename));
nii.hdr.dime.dim(5) = length(flist);
if nii.hdr.dime.dim(1) < 4
nii.hdr.dime.dim(1) = 4;
end
hdr = nii.hdr;
filetype = nii.filetype;
if isfield(nii,'ext') & ~isempty(nii.ext)
ext = nii.ext;
[ext, esize_total] = verify_nii_ext(ext);
else
ext = [];
end
switch double(hdr.dime.datatype),
case 1,
hdr.dime.bitpix = int16(1 ); precision = 'ubit1';
case 2,
hdr.dime.bitpix = int16(8 ); precision = 'uint8';
case 4,
hdr.dime.bitpix = int16(16); precision = 'int16';
case 8,
hdr.dime.bitpix = int16(32); precision = 'int32';
case 16,
hdr.dime.bitpix = int16(32); precision = 'float32';
case 32,
hdr.dime.bitpix = int16(64); precision = 'float32';
case 64,
hdr.dime.bitpix = int16(64); precision = 'float64';
case 128,
hdr.dime.bitpix = int16(24); precision = 'uint8';
case 256
hdr.dime.bitpix = int16(8 ); precision = 'int8';
case 512
hdr.dime.bitpix = int16(16); precision = 'uint16';
case 768
hdr.dime.bitpix = int16(32); precision = 'uint32';
case 1024
hdr.dime.bitpix = int16(64); precision = 'int64';
case 1280
hdr.dime.bitpix = int16(64); precision = 'uint64';
case 1792,
hdr.dime.bitpix = int16(128); precision = 'float64';
otherwise
error('This datatype is not supported');
end
if filetype == 2
fid = fopen(sprintf('%s.nii',fileprefix),'w');
if fid < 0,
msg = sprintf('Cannot open file %s.nii.',fileprefix);
error(msg);
end
hdr.dime.vox_offset = 352;
if ~isempty(ext)
hdr.dime.vox_offset = hdr.dime.vox_offset + esize_total;
end
hdr.hist.magic = 'n+1';
save_untouch_nii_hdr(hdr, fid);
if ~isempty(ext)
save_nii_ext(ext, fid);
end
elseif filetype == 1
fid = fopen(sprintf('%s.hdr',fileprefix),'w');
if fid < 0,
msg = sprintf('Cannot open file %s.hdr.',fileprefix);
error(msg);
end
hdr.dime.vox_offset = 0;
hdr.hist.magic = 'ni1';
save_untouch_nii_hdr(hdr, fid);
if ~isempty(ext)
save_nii_ext(ext, fid);
end
fclose(fid);
fid = fopen(sprintf('%s.img',fileprefix),'w');
else
fid = fopen(sprintf('%s.hdr',fileprefix),'w');
if fid < 0,
msg = sprintf('Cannot open file %s.hdr.',fileprefix);
error(msg);
end
save_untouch0_nii_hdr(hdr, fid);
fclose(fid);
fid = fopen(sprintf('%s.img',fileprefix),'w');
end
if filetype == 2 & isempty(ext)
skip_bytes = double(hdr.dime.vox_offset) - 348;
else
skip_bytes = 0;
end
if skip_bytes
fwrite(fid, zeros(1,skip_bytes), 'uint8');
end
glmax = -inf;
glmin = inf;
for i = 1:length(flist)
nii = load_untouch_nii(fullfile(scan_path,flist{i}));
if double(hdr.dime.datatype) == 128
% RGB planes are expected to be in the 4th dimension of nii.img
%
if(size(nii.img,4)~=3)
error(['The NII structure does not appear to have 3 RGB color planes in the 4th dimension']);
end
nii.img = permute(nii.img, [4 1 2 3 5 6 7 8]);
end
% For complex float32 or complex float64, voxel values
% include [real, imag]
%
if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
real_img = real(nii.img(:))';
nii.img = imag(nii.img(:))';
nii.img = [real_img; nii.img];
end
if nii.hdr.dime.glmax > glmax
glmax = nii.hdr.dime.glmax;
end
if nii.hdr.dime.glmin < glmin
glmin = nii.hdr.dime.glmin;
end
fwrite(fid, nii.img, precision);
end
hdr.dime.glmax = round(glmax);
hdr.dime.glmin = round(glmin);
if filetype == 2
fseek(fid, 140, 'bof');
fwrite(fid, hdr.dime.glmax, 'int32');
fwrite(fid, hdr.dime.glmin, 'int32');
elseif filetype == 1
fid2 = fopen(sprintf('%s.hdr',fileprefix),'w');
if fid2 < 0,
msg = sprintf('Cannot open file %s.hdr.',fileprefix);
error(msg);
end
save_untouch_nii_hdr(hdr, fid2);
if ~isempty(ext)
save_nii_ext(ext, fid2);
end
fclose(fid2);
else
fid2 = fopen(sprintf('%s.hdr',fileprefix),'w');
if fid2 < 0,
msg = sprintf('Cannot open file %s.hdr.',fileprefix);
error(msg);
end
save_untouch0_nii_hdr(hdr, fid2);
fclose(fid2);
end
fclose(fid);
% gzip output file if requested
%
if exist('gzFile', 'var')
if filetype == 1
gzip([fileprefix, '.img']);
delete([fileprefix, '.img']);
gzip([fileprefix, '.hdr']);
delete([fileprefix, '.hdr']);
elseif filetype == 2
gzip([fileprefix, '.nii']);
delete([fileprefix, '.nii']);
end;
end;
return; % collapse_nii_scan
NIfTI_20140122/examples.txt 0 → 100644
View file @ ce12fa61
- Examples to load, make and save a nii struct:
To load Analyze data or NIFTI data to a structure:
nii = load_nii(NIFTI_file_name, [img_idx], [old_RGB24]);
img_idx is a numerical array of image indices along the temporal
axis, which is only available in NIFTI data. After you specify
img_idx, only those images indexed by img_idx will be loaded. If
there is no img_idx or img_idx is empty, all available images
will be loaded.
For RGB image, most people use RGB triple sequentially for each
voxel, like [R1 G1 B1 R2 G2 B2 ...]. However, some program like
Analyze 6.0 developed by AnalyzeDirect uses old RGB24, in a way
like [R1 R2 ... G1 G2 ... B1 B2 ...] for each slices. In this
case, you can set old_RGB24 flag to 1 and load data correctly:
nii = load_nii(NIFTI_file_name, [], 1);
To get a total number of images along the temporal axis:
num_scan = get_nii_frame(NIFTI_file_name);
You can also load the header extension if it exists:
nii.ext = load_nii_ext(NIFTI_file_name);
You can just load the Analyze or NIFTI header:
(header contains: hk, dime, and hist)
hdr = load_nii_hdr(NIFTI_file_name);
You can also save the structure to a new file:
(header extension will be saved if there is nii.ext structure)
save_nii(nii, NIFTI_file_name);
To make the structure from any 3D (or 4D) data:
img = rand(91,109,91); or
img = rand(64,64,21,18);
nii = make_nii(img [, voxel_size, origin, datatype] );
Use "help load_nii", "help save_nii", "help make_nii" etc.
to get more detail information.
- Examples to plot a nii struct:
(More detail descriptions are available on top of "view_nii.m")
Simple way to plot a nii struct:
view_nii(nii);
The default colormap will use the Gray if all data values are
non-negative; otherwise, the default colormap will use BiPolar.
You can choose other colormap, including customized colormap
from panel.
To imbed the plot into your existing figure:
h = gcf;
opt.command = 'init';
opt.setarea = [0.3 0.1 0.6 0.8];
view_nii(h, nii, opt);
To add a colorbar:
opt.usecolorbar = 1;
view_nii(gcf, opt);
Here, opt.command is implicitly set to 'update'.
To display in real aspect ratio:
opt.usestretch = 0;
view_nii(gcf, opt);
If you want the data value to be directly used as the index
of colormap, instead of scale to the whole colormap:
opt.useimagesc = 0;
view_nii(gcf, opt);
If you modified the data value without changing the dimension,
voxel_size, and origin, you can update the display by:
opt.command = 'updateimg';
view_nii(gcf, nii.img, opt);
If the data is completely different, display can be updated by:
opt.command = 'updatenii';
view_nii(gcf, nii, opt);
This is an example to plot EEG source imaging on top of T1 background:
1. download overlay.zip and unzip it from:
http://www.rotman-baycrest.on.ca/~jimmy/NIFTI/overlay.zip
2. T1 = load_nii('T1.nii');
3. EEG = load_nii('EEG.nii');
4. option.setvalue.idx = find(EEG.img);
5. option.setvalue.val = EEG.img(option.setvalue.idx);
6. option.useinterp = 1;
7. option.setviewpoint = [62 48 46];
8. view_nii(T1, option);
- Contrast and Brightness are available under Gray and Bipolar colormap:
Increase contrast in Gray colormap will make high end values
more distinguishable by sacrificing the low end values; The
minimum contrast (default) will display the whole range.
Increase or decrease contrast in BiPolar colormap will shift
the distinguishable position for both positive and negative
values.
Increase or decrease brightness in Gray colormap will shift
the distinguishable position.
Increase or decrease brightness in BiPolar colormap will make
both positive and negative values more distinguishable.
- Required files:
All files in this package.
NIfTI_20140122/expand_nii_scan.m 0 → 100644
View file @ ce12fa61
% Expand a multiple-scan NIFTI file into multiple single-scan NIFTI files
%
% Usage: expand_nii_scan(multi_scan_filename, [img_idx], [path_to_save])
%
% NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function expand_nii_scan(filename, img_idx, newpath)
v = version;
% Check file extension. If .gz, unpack it into temp folder
%
if length(filename) > 2 & strcmp(filename(end-2:end), '.gz')
if ~strcmp(filename(end-6:end), '.img.gz') & ...
~strcmp(filename(end-6:end), '.hdr.gz') & ...
~strcmp(filename(end-6:end), '.nii.gz')
error('Please check filename.');
end
if str2num(v(1:3)) < 7.1 | ~usejava('jvm')
error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.');
else
gzFile = 1;
end
end
if ~exist('newpath','var') | isempty(newpath), newpath = pwd; end
if ~exist('img_idx','var') | isempty(img_idx), img_idx = 1:get_nii_frame(filename); end
for i=img_idx
nii_i = load_untouch_nii(filename, i);
fn = [nii_i.fileprefix '_' sprintf('%04d',i)];
pnfn = fullfile(newpath, fn);
if exist('gzFile', 'var')
pnfn = [pnfn '.nii.gz'];
end
save_untouch_nii(nii_i, pnfn);
end
return; % expand_nii_scan
NIfTI_20140122/extra_nii_hdr.m 0 → 100644
View file @ ce12fa61
% Decode extra NIFTI header information into hdr.extra
%
% Usage: hdr = extra_nii_hdr(hdr)
%
% hdr can be obtained from load_nii_hdr
%
% NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function hdr = extra_nii_hdr(hdr)
switch hdr.dime.datatype
case 1
extra.NIFTI_DATATYPES = 'DT_BINARY';
case 2
extra.NIFTI_DATATYPES = 'DT_UINT8';
case 4
extra.NIFTI_DATATYPES = 'DT_INT16';
case 8
extra.NIFTI_DATATYPES = 'DT_INT32';
case 16
extra.NIFTI_DATATYPES = 'DT_FLOAT32';
case 32
extra.NIFTI_DATATYPES = 'DT_COMPLEX64';
case 64
extra.NIFTI_DATATYPES = 'DT_FLOAT64';
case 128
extra.NIFTI_DATATYPES = 'DT_RGB24';
case 256
extra.NIFTI_DATATYPES = 'DT_INT8';
case 512
extra.NIFTI_DATATYPES = 'DT_UINT16';
case 768
extra.NIFTI_DATATYPES = 'DT_UINT32';
case 1024
extra.NIFTI_DATATYPES = 'DT_INT64';
case 1280
extra.NIFTI_DATATYPES = 'DT_UINT64';
case 1536
extra.NIFTI_DATATYPES = 'DT_FLOAT128';
case 1792
extra.NIFTI_DATATYPES = 'DT_COMPLEX128';
case 2048
extra.NIFTI_DATATYPES = 'DT_COMPLEX256';
otherwise
extra.NIFTI_DATATYPES = 'DT_UNKNOWN';
end
switch hdr.dime.intent_code
case 2
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CORREL';
case 3
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TTEST';
case 4
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_FTEST';
case 5
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_ZSCORE';
case 6
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHISQ';
case 7
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_BETA';
case 8
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_BINOM';
case 9
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_GAMMA';
case 10
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_POISSON';
case 11
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NORMAL';
case 12
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_FTEST_NONC';
case 13
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHISQ_NONC';
case 14
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOGISTIC';
case 15
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LAPLACE';
case 16
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_UNIFORM';
case 17
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TTEST_NONC';
case 18
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_WEIBULL';
case 19
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_CHI';
case 20
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_INVGAUSS';
case 21
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_EXTVAL';
case 22
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_PVAL';
case 23
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOGPVAL';
case 24
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LOG10PVAL';
case 1001
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_ESTIMATE';
case 1002
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_LABEL';
case 1003
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NEURONAME';
case 1004
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_GENMATRIX';
case 1005
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_SYMMATRIX';
case 1006
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_DISPVECT';
case 1007
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_VECTOR';
case 1008
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_POINTSET';
case 1009
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_TRIANGLE';
case 1010
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_QUATERNION';
case 1011
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_DIMLESS';
otherwise
extra.NIFTI_INTENT_CODES = 'NIFTI_INTENT_NONE';
end
extra.NIFTI_INTENT_NAMES = hdr.hist.intent_name;
if hdr.hist.sform_code > 0
switch hdr.hist.sform_code
case 1
extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_SCANNER_ANAT';
case 2
extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_ALIGNED_ANAT';
case 3
extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_TALAIRACH';
case 4
extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_MNI_152';
otherwise
extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
end
extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
elseif hdr.hist.qform_code > 0
extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
switch hdr.hist.qform_code
case 1
extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_SCANNER_ANAT';
case 2
extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_ALIGNED_ANAT';
case 3
extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_TALAIRACH';
case 4
extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_MNI_152';
otherwise
extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
end
else
extra.NIFTI_SFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
extra.NIFTI_QFORM_CODES = 'NIFTI_XFORM_UNKNOWN';
end
switch bitand(hdr.dime.xyzt_units, 7) % mask with 0x07
case 1
extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_METER';
case 2
extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_MM'; % millimeter
case 3
extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_MICRO';
otherwise
extra.NIFTI_SPACE_UNIT = 'NIFTI_UNITS_UNKNOWN';
end
switch bitand(hdr.dime.xyzt_units, 56) % mask with 0x38
case 8
extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_SEC';
case 16
extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_MSEC';
case 24
extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_USEC'; % microsecond
otherwise
extra.NIFTI_TIME_UNIT = 'NIFTI_UNITS_UNKNOWN';
end
switch hdr.dime.xyzt_units
case 32
extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_HZ';
case 40
extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_PPM'; % part per million
case 48
extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_RADS'; % radians per second
otherwise
extra.NIFTI_SPECTRAL_UNIT = 'NIFTI_UNITS_UNKNOWN';
end
% MRI-specific spatial and temporal information
%
dim_info = hdr.hk.dim_info;
extra.NIFTI_FREQ_DIM = bitand(dim_info, 3);
extra.NIFTI_PHASE_DIM = bitand(bitshift(dim_info, -2), 3);
extra.NIFTI_SLICE_DIM = bitand(bitshift(dim_info, -4), 3);
% Check slice code
%
switch hdr.dime.slice_code
case 1
extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_SEQ_INC'; % sequential increasing
case 2
extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_SEQ_DEC'; % sequential decreasing
case 3
extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_INC'; % alternating increasing
case 4
extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_DEC'; % alternating decreasing
case 5
extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_INC2'; % ALT_INC # 2
case 6
extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_ALT_DEC2'; % ALT_DEC # 2
otherwise
extra.NIFTI_SLICE_ORDER = 'NIFTI_SLICE_UNKNOWN';
end
% Check NIFTI version
%
if ~isempty(hdr.hist.magic) & strcmp(hdr.hist.magic(1),'n') & ...
( strcmp(hdr.hist.magic(2),'i') | strcmp(hdr.hist.magic(2),'+') ) & ...
str2num(hdr.hist.magic(3)) >= 1 & str2num(hdr.hist.magic(3)) <= 9
extra.NIFTI_VERSION = str2num(hdr.hist.magic(3));
else
extra.NIFTI_VERSION = 0;
end
% Check if data stored in the same file (*.nii) or separate
% files (*.hdr/*.img)
%
if isempty(hdr.hist.magic)
extra.NIFTI_ONEFILE = 0;
else
extra.NIFTI_ONEFILE = strcmp(hdr.hist.magic(2), '+');
end
% Swap has been taken care of by checking whether sizeof_hdr is
% 348 (machine is 'ieee-le' or 'ieee-be' etc)
%
% extra.NIFTI_NEEDS_SWAP = (hdr.dime.dim(1) < 0 | hdr.dime.dim(1) > 7);
% Check NIFTI header struct contains a 5th (vector) dimension
%
if hdr.dime.dim(1) > 4 & hdr.dime.dim(6) > 1
extra.NIFTI_5TH_DIM = hdr.dime.dim(6);
else
extra.NIFTI_5TH_DIM = 0;
end
hdr.extra = extra;
return; % extra_nii_hdr
NIfTI_20140122/flip_lr.m 0 → 100644
View file @ ce12fa61
% When you load any ANALYZE or NIfTI file with 'load_nii.m', and view
% it with 'view_nii.m', you may find that the image is L-R flipped.
% This is because of the confusion of radiological and neurological
% convention in the medical image before NIfTI format is adopted. You
% can find more details from:
%
% http://www.rotman-baycrest.on.ca/~jimmy/UseANALYZE.htm
%
% Sometime, people even want to convert RAS (standard orientation) back
% to LAS orientation to satisfy the legend programs or processes. This
% program is only written for those purpose. So PLEASE BE VERY CAUTIOUS
% WHEN USING THIS 'FLIP_LR.M' PROGRAM.
%
% With 'flip_lr.m', you can convert any ANALYZE or NIfTI (no matter
% 3D or 4D) file to a flipped NIfTI file. This is implemented simply
% by flipping the affine matrix in the NIfTI header. Since the L-R
% orientation is determined there, so the image will be flipped.
%
% Usage: flip_lr(original_fn, flipped_fn, [old_RGB],[tolerance],[preferredForm])
%
% original_fn - filename of the original ANALYZE or NIfTI (3D or 4D) file
%
% flipped_fn - filename of the L-R flipped NIfTI file
%
% old_RGB (optional) - a scale number to tell difference of new RGB24
% from old RGB24. New RGB24 uses RGB triple sequentially for each
% voxel, like [R1 G1 B1 R2 G2 B2 ...]. Analyze 6.0 from AnalyzeDirect
% uses old RGB24, in a way like [R1 R2 ... G1 G2 ... B1 B2 ...] for
% each slices. If the image that you view is garbled, try to set
% old_RGB variable to 1 and try again, because it could be in
% old RGB24. It will be set to 0, if it is default or empty.
%
% tolerance (optional) - distortion allowed for non-orthogonal rotation
% or shearing in NIfTI affine matrix. It will be set to 0.1 (10%),
% if it is default or empty.
%
% preferredForm (optional) - selects which transformation from voxels
% to RAS coordinates; values are s,q,S,Q. Lower case s,q indicate
% "prefer sform or qform, but use others if preferred not present".
% Upper case indicate the program is forced to use the specificied
% tranform or fail loading. 'preferredForm' will be 's', if it is
% default or empty. - Jeff Gunter
%
% Example: flip_lr('avg152T1_LR_nifti.nii', 'flipped_lr.nii');
% flip_lr('avg152T1_RL_nifti.nii', 'flipped_rl.nii');
%
% You will find that 'avg152T1_LR_nifti.nii' and 'avg152T1_RL_nifti.nii'
% are the same, and 'flipped_lr.nii' and 'flipped_rl.nii' are also the
% the same, but they are L-R flipped from 'avg152T1_*'.
%
% NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function flip_lr(original_fn, flipped_fn, old_RGB, tolerance, preferredForm)
if ~exist('original_fn','var') | ~exist('flipped_fn','var')
error('Usage: flip_lr(original_fn, flipped_fn, [old_RGB],[tolerance])');
end
if ~exist('old_RGB','var') | isempty(old_RGB)
old_RGB = 0;
end
if ~exist('tolerance','var') | isempty(tolerance)
tolerance = 0.1;
end
if ~exist('preferredForm','var') | isempty(preferredForm)
preferredForm= 's'; % Jeff
end
nii = load_nii(original_fn, [], [], [], [], old_RGB, tolerance, preferredForm);
M = diag(nii.hdr.dime.pixdim(2:5));
M(1:3,4) = -M(1:3,1:3)*(nii.hdr.hist.originator(1:3)-1)';
M(1,:) = -1*M(1,:);
nii.hdr.hist.sform_code = 1;
nii.hdr.hist.srow_x = M(1,:);
nii.hdr.hist.srow_y = M(2,:);
nii.hdr.hist.srow_z = M(3,:);
save_nii(nii, flipped_fn);
return; % flip_lr
NIfTI_20140122/get_nii_frame.m 0 → 100644
View file @ ce12fa61
% Return time frame of a NIFTI dataset. Support both *.nii and
% *.hdr/*.img file extension. If file extension is not provided,
% *.hdr/*.img will be used as default.
%
% It is a lightweighted "load_nii_hdr", and is equivalent to
% hdr.dime.dim(5)
%
% Usage: [ total_scan ] = get_nii_frame(filename)
%
% filename - NIFTI file name.
%
% Returned values:
%
% total_scan - total number of image scans for the time frame
%
% NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function [ total_scan ] = get_nii_frame(filename)
if ~exist('filename','var'),
error('Usage: [ total_scan ] = get_nii_frame(filename)');
end
v = version;
% Check file extension. If .gz, unpack it into temp folder
%
if length(filename) > 2 & strcmp(filename(end-2:end), '.gz')
if ~strcmp(filename(end-6:end), '.img.gz') & ...
~strcmp(filename(end-6:end), '.hdr.gz') & ...
~strcmp(filename(end-6:end), '.nii.gz')
error('Please check filename.');
end
if str2num(v(1:3)) < 7.1 | ~usejava('jvm')
error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.');
elseif strcmp(filename(end-6:end), '.img.gz')
filename1 = filename;
filename2 = filename;
filename2(end-6:end) = '';
filename2 = [filename2, '.hdr.gz'];
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename1 = gunzip(filename1, tmpDir);
filename2 = gunzip(filename2, tmpDir);
filename = char(filename1); % convert from cell to string
elseif strcmp(filename(end-6:end), '.hdr.gz')
filename1 = filename;
filename2 = filename;
filename2(end-6:end) = '';
filename2 = [filename2, '.img.gz'];
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename1 = gunzip(filename1, tmpDir);
filename2 = gunzip(filename2, tmpDir);
filename = char(filename1); % convert from cell to string
elseif strcmp(filename(end-6:end), '.nii.gz')
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename = gunzip(filename, tmpDir);
filename = char(filename); % convert from cell to string
end
end
fileprefix = filename;
machine = 'ieee-le';
new_ext = 0;
if findstr('.nii',fileprefix) & strcmp(fileprefix(end-3:end), '.nii')
new_ext = 1;
fileprefix(end-3:end)='';
end
if findstr('.hdr',fileprefix) & strcmp(fileprefix(end-3:end), '.hdr')
fileprefix(end-3:end)='';
end
if findstr('.img',fileprefix) & strcmp(fileprefix(end-3:end), '.img')
fileprefix(end-3:end)='';
end
if new_ext
fn = sprintf('%s.nii',fileprefix);
if ~exist(fn)
msg = sprintf('Cannot find file "%s.nii".', fileprefix);
error(msg);
end
else
fn = sprintf('%s.hdr',fileprefix);
if ~exist(fn)
msg = sprintf('Cannot find file "%s.hdr".', fileprefix);
error(msg);
end
end
fid = fopen(fn,'r',machine);
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
else
hdr = read_header(fid);
fclose(fid);
end
if hdr.sizeof_hdr ~= 348
% first try reading the opposite endian to 'machine'
switch machine,
case 'ieee-le', machine = 'ieee-be';
case 'ieee-be', machine = 'ieee-le';
end
fid = fopen(fn,'r',machine);
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
else
hdr = read_header(fid);
fclose(fid);
end
end
if hdr.sizeof_hdr ~= 348
% Now throw an error
msg = sprintf('File "%s" is corrupted.',fn);
error(msg);
end
total_scan = hdr.dim(5);
% Clean up after gunzip
%
if exist('gzFileName', 'var')
rmdir(tmpDir,'s');
end
return; % get_nii_frame
%---------------------------------------------------------------------
function [ dsr ] = read_header(fid)
fseek(fid,0,'bof');
dsr.sizeof_hdr = fread(fid,1,'int32')'; % should be 348!
fseek(fid,40,'bof');
dsr.dim = fread(fid,8,'int16')';
return; % read_header
NIfTI_20140122/license.txt 0 → 100644
View file @ ce12fa61
Copyright (c) 2014, Jimmy Shen
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
NIfTI_20140122/load_nii.m 0 → 100644
View file @ ce12fa61
% Load NIFTI or ANALYZE dataset. Support both *.nii and *.hdr/*.img
% file extension. If file extension is not provided, *.hdr/*.img will
% be used as default.
%
% A subset of NIFTI transform is included. For non-orthogonal rotation,
% shearing etc., please use 'reslice_nii.m' to reslice the NIFTI file.
% It will not cause negative effect, as long as you remember not to do
% slice time correction after reslicing the NIFTI file. Output variable
% nii will be in RAS orientation, i.e. X axis from Left to Right,
% Y axis from Posterior to Anterior, and Z axis from Inferior to
% Superior.
%
% Usage: nii = load_nii(filename, [img_idx], [dim5_idx], [dim6_idx], ...
% [dim7_idx], [old_RGB], [tolerance], [preferredForm])
%
% filename - NIFTI or ANALYZE file name.
%
% img_idx (optional) - a numerical array of 4th dimension indices,
% which is the indices of image scan volume. The number of images
% scan volumes can be obtained from get_nii_frame.m, or simply
% hdr.dime.dim(5). Only the specified volumes will be loaded.
% All available image volumes will be loaded, if it is default or
% empty.
%
% dim5_idx (optional) - a numerical array of 5th dimension indices.
% Only the specified range will be loaded. All available range
% will be loaded, if it is default or empty.
%
% dim6_idx (optional) - a numerical array of 6th dimension indices.
% Only the specified range will be loaded. All available range
% will be loaded, if it is default or empty.
%
% dim7_idx (optional) - a numerical array of 7th dimension indices.
% Only the specified range will be loaded. All available range
% will be loaded, if it is default or empty.
%
% old_RGB (optional) - a scale number to tell difference of new RGB24
% from old RGB24. New RGB24 uses RGB triple sequentially for each
% voxel, like [R1 G1 B1 R2 G2 B2 ...]. Analyze 6.0 from AnalyzeDirect
% uses old RGB24, in a way like [R1 R2 ... G1 G2 ... B1 B2 ...] for
% each slices. If the image that you view is garbled, try to set
% old_RGB variable to 1 and try again, because it could be in
% old RGB24. It will be set to 0, if it is default or empty.
%
% tolerance (optional) - distortion allowed in the loaded image for any
% non-orthogonal rotation or shearing of NIfTI affine matrix. If
% you set 'tolerance' to 0, it means that you do not allow any
% distortion. If you set 'tolerance' to 1, it means that you do
% not care any distortion. The image will fail to be loaded if it
% can not be tolerated. The tolerance will be set to 0.1 (10%), if
% it is default or empty.
%
% preferredForm (optional) - selects which transformation from voxels
% to RAS coordinates; values are s,q,S,Q. Lower case s,q indicate
% "prefer sform or qform, but use others if preferred not present".
% Upper case indicate the program is forced to use the specificied
% tranform or fail loading. 'preferredForm' will be 's', if it is
% default or empty. - Jeff Gunter
%
% Returned values:
%
% nii structure:
%
% hdr - struct with NIFTI header fields.
%
% filetype - Analyze format .hdr/.img (0);
% NIFTI .hdr/.img (1);
% NIFTI .nii (2)
%
% fileprefix - NIFTI filename without extension.
%
% machine - machine string variable.
%
% img - 3D (or 4D) matrix of NIFTI data.
%
% original - the original header before any affine transform.
%
% Part of this file is copied and modified from:
% http://www.mathworks.com/matlabcentral/fileexchange/1878-mri-analyze-tools
%
% NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function nii = load_nii(filename, img_idx, dim5_idx, dim6_idx, dim7_idx, ...
old_RGB, tolerance, preferredForm)
if ~exist('filename','var')
error('Usage: nii = load_nii(filename, [img_idx], [dim5_idx], [dim6_idx], [dim7_idx], [old_RGB], [tolerance], [preferredForm])');
end
if ~exist('img_idx','var') | isempty(img_idx)
img_idx = [];
end
if ~exist('dim5_idx','var') | isempty(dim5_idx)
dim5_idx = [];
end
if ~exist('dim6_idx','var') | isempty(dim6_idx)
dim6_idx = [];
end
if ~exist('dim7_idx','var') | isempty(dim7_idx)
dim7_idx = [];
end
if ~exist('old_RGB','var') | isempty(old_RGB)
old_RGB = 0;
end
if ~exist('tolerance','var') | isempty(tolerance)
tolerance = 0.1; % 10 percent
end
if ~exist('preferredForm','var') | isempty(preferredForm)
preferredForm= 's'; % Jeff
end
v = version;
% Check file extension. If .gz, unpack it into temp folder
%
if length(filename) > 2 & strcmp(filename(end-2:end), '.gz')
if ~strcmp(filename(end-6:end), '.img.gz') & ...
~strcmp(filename(end-6:end), '.hdr.gz') & ...
~strcmp(filename(end-6:end), '.nii.gz')
error('Please check filename.');
end
if str2num(v(1:3)) < 7.1 | ~usejava('jvm')
error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.');
elseif strcmp(filename(end-6:end), '.img.gz')
filename1 = filename;
filename2 = filename;
filename2(end-6:end) = '';
filename2 = [filename2, '.hdr.gz'];
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename1 = gunzip(filename1, tmpDir);
filename2 = gunzip(filename2, tmpDir);
filename = char(filename1); % convert from cell to string
elseif strcmp(filename(end-6:end), '.hdr.gz')
filename1 = filename;
filename2 = filename;
filename2(end-6:end) = '';
filename2 = [filename2, '.img.gz'];
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename1 = gunzip(filename1, tmpDir);
filename2 = gunzip(filename2, tmpDir);
filename = char(filename1); % convert from cell to string
elseif strcmp(filename(end-6:end), '.nii.gz')
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename = gunzip(filename, tmpDir);
filename = char(filename); % convert from cell to string
end
end
% Read the dataset header
%
[nii.hdr,nii.filetype,nii.fileprefix,nii.machine] = load_nii_hdr(filename);
% Read the header extension
%
% nii.ext = load_nii_ext(filename);
% Read the dataset body
%
[nii.img,nii.hdr] = load_nii_img(nii.hdr,nii.filetype,nii.fileprefix, ...
nii.machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB);
% Perform some of sform/qform transform
%
nii = xform_nii(nii, tolerance, preferredForm);
% Clean up after gunzip
%
if exist('gzFileName', 'var')
% fix fileprefix so it doesn't point to temp location
%
nii.fileprefix = gzFileName(1:end-7);
rmdir(tmpDir,'s');
end
return % load_nii
NIfTI_20140122/load_nii_ext.m 0 → 100644
View file @ ce12fa61
% Load NIFTI header extension after its header is loaded using load_nii_hdr.
%
% Usage: ext = load_nii_ext(filename)
%
% filename - NIFTI file name.
%
% Returned values:
%
% ext - Structure of NIFTI header extension, which includes num_ext,
% and all the extended header sections in the header extension.
% Each extended header section will have its esize, ecode, and
% edata, where edata can be plain text, xml, or any raw data
% that was saved in the extended header section.
%
% NIFTI data format can be found on: http://nifti.nimh.nih.gov
%
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function ext = load_nii_ext(filename)
if ~exist('filename','var'),
error('Usage: ext = load_nii_ext(filename)');
end
v = version;
% Check file extension. If .gz, unpack it into temp folder
%
if length(filename) > 2 & strcmp(filename(end-2:end), '.gz')
if ~strcmp(filename(end-6:end), '.img.gz') & ...
~strcmp(filename(end-6:end), '.hdr.gz') & ...
~strcmp(filename(end-6:end), '.nii.gz')
error('Please check filename.');
end
if str2num(v(1:3)) < 7.1 | ~usejava('jvm')
error('Please use MATLAB 7.1 (with java) and above, or run gunzip outside MATLAB.');
elseif strcmp(filename(end-6:end), '.img.gz')
filename1 = filename;
filename2 = filename;
filename2(end-6:end) = '';
filename2 = [filename2, '.hdr.gz'];
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename1 = gunzip(filename1, tmpDir);
filename2 = gunzip(filename2, tmpDir);
filename = char(filename1); % convert from cell to string
elseif strcmp(filename(end-6:end), '.hdr.gz')
filename1 = filename;
filename2 = filename;
filename2(end-6:end) = '';
filename2 = [filename2, '.img.gz'];
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename1 = gunzip(filename1, tmpDir);
filename2 = gunzip(filename2, tmpDir);
filename = char(filename1); % convert from cell to string
elseif strcmp(filename(end-6:end), '.nii.gz')
tmpDir = tempname;
mkdir(tmpDir);
gzFileName = filename;
filename = gunzip(filename, tmpDir);
filename = char(filename); % convert from cell to string
end
end
machine = 'ieee-le';
new_ext = 0;
if findstr('.nii',filename) & strcmp(filename(end-3:end), '.nii')
new_ext = 1;
filename(end-3:end)='';
end
if findstr('.hdr',filename) & strcmp(filename(end-3:end), '.hdr')
filename(end-3:end)='';
end
if findstr('.img',filename) & strcmp(filename(end-3:end), '.img')
filename(end-3:end)='';
end
if new_ext
fn = sprintf('%s.nii',filename);
if ~exist(fn)
msg = sprintf('Cannot find file "%s.nii".', filename);
error(msg);
end
else
fn = sprintf('%s.hdr',filename);
if ~exist(fn)
msg = sprintf('Cannot find file "%s.hdr".', filename);
error(msg);
end
end
fid = fopen(fn,'r',machine);
vox_offset = 0;
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
else
fseek(fid,0,'bof');
if fread(fid,1,'int32') == 348
if new_ext
fseek(fid,108,'bof');
vox_offset = fread(fid,1,'float32');
end
ext = read_extension(fid, vox_offset);
fclose(fid);
else
fclose(fid);
% first try reading the opposite endian to 'machine'
%
switch machine,
case 'ieee-le', machine = 'ieee-be';
case 'ieee-be', machine = 'ieee-le';
end
fid = fopen(fn,'r',machine);
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
else
fseek(fid,0,'bof');
if fread(fid,1,'int32') ~= 348
% Now throw an error
%
msg = sprintf('File "%s" is corrupted.',fn);
error(msg);
end
if new_ext
fseek(fid,108,'bof');
vox_offset = fread(fid,1,'float32');
end
ext = read_extension(fid, vox_offset);
fclose(fid);
end
end
end
% Clean up after gunzip
%
if exist('gzFileName', 'var')
rmdir(tmpDir,'s');
end
return % load_nii_ext
%---------------------------------------------------------------------
function ext = read_extension(fid, vox_offset)
ext = [];
if vox_offset
end_of_ext = vox_offset;
else
fseek(fid, 0, 'eof');
end_of_ext = ftell(fid);
end
if end_of_ext > 352
fseek(fid, 348, 'bof');
ext.extension = fread(fid,4)';
end
if isempty(ext) | ext.extension(1) == 0
ext = [];
return;
end
i = 1;
while(ftell(fid) < end_of_ext)
ext.section(i).esize = fread(fid,1,'int32');
ext.section(i).ecode = fread(fid,1,'int32');
ext.section(i).edata = char(fread(fid,ext.section(i).esize-8)');
i = i + 1;
end
ext.num_ext = length(ext.section);
return % read_extension
NIfTI_20140122/load_nii_hdr.m 0 → 100644
View file @ ce12fa61
% internal function
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
function [hdr, filetype, fileprefix, machine] = load_nii_hdr(fileprefix)
if ~exist('fileprefix','var'),
error('Usage: [hdr, filetype, fileprefix, machine] = load_nii_hdr(filename)');
end
machine = 'ieee-le';
new_ext = 0;
if findstr('.nii',fileprefix) & strcmp(fileprefix(end-3:end), '.nii')
new_ext = 1;
fileprefix(end-3:end)='';
end
if findstr('.hdr',fileprefix) & strcmp(fileprefix(end-3:end), '.hdr')
fileprefix(end-3:end)='';
end
if findstr('.img',fileprefix) & strcmp(fileprefix(end-3:end), '.img')
fileprefix(end-3:end)='';
end
if new_ext
fn = sprintf('%s.nii',fileprefix);
if ~exist(fn)
msg = sprintf('Cannot find file "%s.nii".', fileprefix);
error(msg);
end
else
fn = sprintf('%s.hdr',fileprefix);
if ~exist(fn)
msg = sprintf('Cannot find file "%s.hdr".', fileprefix);
error(msg);
end
end
fid = fopen(fn,'r',machine);
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
else
fseek(fid,0,'bof');
if fread(fid,1,'int32') == 348
hdr = read_header(fid);
fclose(fid);
else
fclose(fid);
% first try reading the opposite endian to 'machine'
%
switch machine,
case 'ieee-le', machine = 'ieee-be';
case 'ieee-be', machine = 'ieee-le';
end
fid = fopen(fn,'r',machine);
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
else
fseek(fid,0,'bof');
if fread(fid,1,'int32') ~= 348
% Now throw an error
%
msg = sprintf('File "%s" is corrupted.',fn);
error(msg);
end
hdr = read_header(fid);
fclose(fid);
end
end
end
if strcmp(hdr.hist.magic, 'n+1')
filetype = 2;
elseif strcmp(hdr.hist.magic, 'ni1')
filetype = 1;
else
filetype = 0;
end
return % load_nii_hdr
%---------------------------------------------------------------------
function [ dsr ] = read_header(fid)
% Original header structures
% struct dsr
% {
% struct header_key hk; /* 0 + 40 */
% struct image_dimension dime; /* 40 + 108 */
% struct data_history hist; /* 148 + 200 */
% }; /* total= 348 bytes*/
dsr.hk = header_key(fid);
dsr.dime = image_dimension(fid);
dsr.hist = data_history(fid);
% For Analyze data format
%
if ~strcmp(dsr.hist.magic, 'n+1') & ~strcmp(dsr.hist.magic, 'ni1')
dsr.hist.qform_code = 0;
dsr.hist.sform_code = 0;
end
return % read_header
%---------------------------------------------------------------------
function [ hk ] = header_key(fid)
fseek(fid,0,'bof');
% Original header structures
% struct header_key /* header key */
% { /* off + size */
% int sizeof_hdr /* 0 + 4 */
% char data_type[10]; /* 4 + 10 */
% char db_name[18]; /* 14 + 18 */
% int extents; /* 32 + 4 */
% short int session_error; /* 36 + 2 */
% char regular; /* 38 + 1 */
% char dim_info; % char hkey_un0; /* 39 + 1 */
% }; /* total=40 bytes */
%
% int sizeof_header Should be 348.
% char regular Must be 'r' to indicate that all images and
% volumes are the same size.
v6 = version;
if str2num(v6(1))<6
directchar = '*char';
else
directchar = 'uchar=>char';
end
hk.sizeof_hdr = fread(fid, 1,'int32')'; % should be 348!
hk.data_type = deblank(fread(fid,10,directchar)');
hk.db_name = deblank(fread(fid,18,directchar)');
hk.extents = fread(fid, 1,'int32')';
hk.session_error = fread(fid, 1,'int16')';
hk.regular = fread(fid, 1,directchar)';
hk.dim_info = fread(fid, 1,'uchar')';
return % header_key
%---------------------------------------------------------------------
function [ dime ] = image_dimension(fid)
% Original header structures
% struct image_dimension
% { /* off + size */
% short int dim[8]; /* 0 + 16 */
% /*
% dim[0] Number of dimensions in database; usually 4.
% dim[1] Image X dimension; number of *pixels* in an image row.
% dim[2] Image Y dimension; number of *pixel rows* in slice.
% dim[3] Volume Z dimension; number of *slices* in a volume.
% dim[4] Time points; number of volumes in database
% */
% float intent_p1; % char vox_units[4]; /* 16 + 4 */
% float intent_p2; % char cal_units[8]; /* 20 + 4 */
% float intent_p3; % char cal_units[8]; /* 24 + 4 */
% short int intent_code; % short int unused1; /* 28 + 2 */
% short int datatype; /* 30 + 2 */
% short int bitpix; /* 32 + 2 */
% short int slice_start; % short int dim_un0; /* 34 + 2 */
% float pixdim[8]; /* 36 + 32 */
% /*
% pixdim[] specifies the voxel dimensions:
% pixdim[1] - voxel width, mm
% pixdim[2] - voxel height, mm
% pixdim[3] - slice thickness, mm
% pixdim[4] - volume timing, in msec
% ..etc
% */
% float vox_offset; /* 68 + 4 */
% float scl_slope; % float roi_scale; /* 72 + 4 */
% float scl_inter; % float funused1; /* 76 + 4 */
% short slice_end; % float funused2; /* 80 + 2 */
% char slice_code; % float funused2; /* 82 + 1 */
% char xyzt_units; % float funused2; /* 83 + 1 */
% float cal_max; /* 84 + 4 */
% float cal_min; /* 88 + 4 */
% float slice_duration; % int compressed; /* 92 + 4 */
% float toffset; % int verified; /* 96 + 4 */
% int glmax; /* 100 + 4 */
% int glmin; /* 104 + 4 */
% }; /* total=108 bytes */
dime.dim = fread(fid,8,'int16')';
dime.intent_p1 = fread(fid,1,'float32')';
dime.intent_p2 = fread(fid,1,'float32')';
dime.intent_p3 = fread(fid,1,'float32')';
dime.intent_code = fread(fid,1,'int16')';
dime.datatype = fread(fid,1,'int16')';
dime.bitpix = fread(fid,1,'int16')';
dime.slice_start = fread(fid,1,'int16')';
dime.pixdim = fread(fid,8,'float32')';
dime.vox_offset = fread(fid,1,'float32')';
dime.scl_slope = fread(fid,1,'float32')';
dime.scl_inter = fread(fid,1,'float32')';
dime.slice_end = fread(fid,1,'int16')';
dime.slice_code = fread(fid,1,'uchar')';
dime.xyzt_units = fread(fid,1,'uchar')';
dime.cal_max = fread(fid,1,'float32')';
dime.cal_min = fread(fid,1,'float32')';
dime.slice_duration = fread(fid,1,'float32')';
dime.toffset = fread(fid,1,'float32')';
dime.glmax = fread(fid,1,'int32')';
dime.glmin = fread(fid,1,'int32')';
return % image_dimension
%---------------------------------------------------------------------
function [ hist ] = data_history(fid)
% Original header structures
% struct data_history
% { /* off + size */
% char descrip[80]; /* 0 + 80 */
% char aux_file[24]; /* 80 + 24 */
% short int qform_code; /* 104 + 2 */
% short int sform_code; /* 106 + 2 */
% float quatern_b; /* 108 + 4 */
% float quatern_c; /* 112 + 4 */
% float quatern_d; /* 116 + 4 */
% float qoffset_x; /* 120 + 4 */
% float qoffset_y; /* 124 + 4 */
% float qoffset_z; /* 128 + 4 */
% float srow_x[4]; /* 132 + 16 */
% float srow_y[4]; /* 148 + 16 */
% float srow_z[4]; /* 164 + 16 */
% char intent_name[16]; /* 180 + 16 */
% char magic[4]; % int smin; /* 196 + 4 */
% }; /* total=200 bytes */
v6 = version;
if str2num(v6(1))<6
directchar = '*char';
else
directchar = 'uchar=>char';
end
hist.descrip = deblank(fread(fid,80,directchar)');
hist.aux_file = deblank(fread(fid,24,directchar)');
hist.qform_code = fread(fid,1,'int16')';
hist.sform_code = fread(fid,1,'int16')';
hist.quatern_b = fread(fid,1,'float32')';
hist.quatern_c = fread(fid,1,'float32')';
hist.quatern_d = fread(fid,1,'float32')';
hist.qoffset_x = fread(fid,1,'float32')';
hist.qoffset_y = fread(fid,1,'float32')';
hist.qoffset_z = fread(fid,1,'float32')';
hist.srow_x = fread(fid,4,'float32')';
hist.srow_y = fread(fid,4,'float32')';
hist.srow_z = fread(fid,4,'float32')';
hist.intent_name = deblank(fread(fid,16,directchar)');
hist.magic = deblank(fread(fid,4,directchar)');
fseek(fid,253,'bof');
hist.originator = fread(fid, 5,'int16')';
return % data_history
NIfTI_20140122/load_nii_img.m 0 → 100644
View file @ ce12fa61
% internal function
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
function [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)
if ~exist('hdr','var') | ~exist('filetype','var') | ~exist('fileprefix','var') | ~exist('machine','var')
error('Usage: [img,hdr] = load_nii_img(hdr,filetype,fileprefix,machine,[img_idx],[dim5_idx],[dim6_idx],[dim7_idx],[old_RGB]);');
end
if ~exist('img_idx','var') | isempty(img_idx) | hdr.dime.dim(5)<1
img_idx = [];
end
if ~exist('dim5_idx','var') | isempty(dim5_idx) | hdr.dime.dim(6)<1
dim5_idx = [];
end
if ~exist('dim6_idx','var') | isempty(dim6_idx) | hdr.dime.dim(7)<1
dim6_idx = [];
end
if ~exist('dim7_idx','var') | isempty(dim7_idx) | hdr.dime.dim(8)<1
dim7_idx = [];
end
if ~exist('old_RGB','var') | isempty(old_RGB)
old_RGB = 0;
end
% check img_idx
%
if ~isempty(img_idx) & ~isnumeric(img_idx)
error('"img_idx" should be a numerical array.');
end
if length(unique(img_idx)) ~= length(img_idx)
error('Duplicate image index in "img_idx"');
end
if ~isempty(img_idx) & (min(img_idx) < 1 | max(img_idx) > hdr.dime.dim(5))
max_range = hdr.dime.dim(5);
if max_range == 1
error(['"img_idx" should be 1.']);
else
range = ['1 ' num2str(max_range)];
error(['"img_idx" should be an integer within the range of [' range '].']);
end
end
% check dim5_idx
%
if ~isempty(dim5_idx) & ~isnumeric(dim5_idx)
error('"dim5_idx" should be a numerical array.');
end
if length(unique(dim5_idx)) ~= length(dim5_idx)
error('Duplicate index in "dim5_idx"');
end
if ~isempty(dim5_idx) & (min(dim5_idx) < 1 | max(dim5_idx) > hdr.dime.dim(6))
max_range = hdr.dime.dim(6);
if max_range == 1
error(['"dim5_idx" should be 1.']);
else
range = ['1 ' num2str(max_range)];
error(['"dim5_idx" should be an integer within the range of [' range '].']);
end
end
% check dim6_idx
%
if ~isempty(dim6_idx) & ~isnumeric(dim6_idx)
error('"dim6_idx" should be a numerical array.');
end
if length(unique(dim6_idx)) ~= length(dim6_idx)
error('Duplicate index in "dim6_idx"');
end
if ~isempty(dim6_idx) & (min(dim6_idx) < 1 | max(dim6_idx) > hdr.dime.dim(7))
max_range = hdr.dime.dim(7);
if max_range == 1
error(['"dim6_idx" should be 1.']);
else
range = ['1 ' num2str(max_range)];
error(['"dim6_idx" should be an integer within the range of [' range '].']);
end
end
% check dim7_idx
%
if ~isempty(dim7_idx) & ~isnumeric(dim7_idx)
error('"dim7_idx" should be a numerical array.');
end
if length(unique(dim7_idx)) ~= length(dim7_idx)
error('Duplicate index in "dim7_idx"');
end
if ~isempty(dim7_idx) & (min(dim7_idx) < 1 | max(dim7_idx) > hdr.dime.dim(8))
max_range = hdr.dime.dim(8);
if max_range == 1
error(['"dim7_idx" should be 1.']);
else
range = ['1 ' num2str(max_range)];
error(['"dim7_idx" should be an integer within the range of [' range '].']);
end
end
[img,hdr] = read_image(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB);
return % load_nii_img
%---------------------------------------------------------------------
function [img,hdr] = read_image(hdr,filetype,fileprefix,machine,img_idx,dim5_idx,dim6_idx,dim7_idx,old_RGB)
switch filetype
case {0, 1}
fn = [fileprefix '.img'];
case 2
fn = [fileprefix '.nii'];
end
fid = fopen(fn,'r',machine);
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
end
% Set bitpix according to datatype
%
% /*Acceptable values for datatype are*/
%
% 0 None (Unknown bit per voxel) % DT_NONE, DT_UNKNOWN
% 1 Binary (ubit1, bitpix=1) % DT_BINARY
% 2 Unsigned char (uchar or uint8, bitpix=8) % DT_UINT8, NIFTI_TYPE_UINT8
% 4 Signed short (int16, bitpix=16) % DT_INT16, NIFTI_TYPE_INT16
% 8 Signed integer (int32, bitpix=32) % DT_INT32, NIFTI_TYPE_INT32
% 16 Floating point (single or float32, bitpix=32) % DT_FLOAT32, NIFTI_TYPE_FLOAT32
% 32 Complex, 2 float32 (Use float32, bitpix=64) % DT_COMPLEX64, NIFTI_TYPE_COMPLEX64
% 64 Double precision (double or float64, bitpix=64) % DT_FLOAT64, NIFTI_TYPE_FLOAT64
% 128 uint8 RGB (Use uint8, bitpix=24) % DT_RGB24, NIFTI_TYPE_RGB24
% 256 Signed char (schar or int8, bitpix=8) % DT_INT8, NIFTI_TYPE_INT8
% 511 Single RGB (Use float32, bitpix=96) % DT_RGB96, NIFTI_TYPE_RGB96
% 512 Unsigned short (uint16, bitpix=16) % DT_UNINT16, NIFTI_TYPE_UNINT16
% 768 Unsigned integer (uint32, bitpix=32) % DT_UNINT32, NIFTI_TYPE_UNINT32
% 1024 Signed long long (int64, bitpix=64) % DT_INT64, NIFTI_TYPE_INT64
% 1280 Unsigned long long (uint64, bitpix=64) % DT_UINT64, NIFTI_TYPE_UINT64
% 1536 Long double, float128 (Unsupported, bitpix=128) % DT_FLOAT128, NIFTI_TYPE_FLOAT128
% 1792 Complex128, 2 float64 (Use float64, bitpix=128) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
% 2048 Complex256, 2 float128 (Unsupported, bitpix=256) % DT_COMPLEX128, NIFTI_TYPE_COMPLEX128
%
switch hdr.dime.datatype
case 1,
hdr.dime.bitpix = 1; precision = 'ubit1';
case 2,
hdr.dime.bitpix = 8; precision = 'uint8';
case 4,
hdr.dime.bitpix = 16; precision = 'int16';
case 8,
hdr.dime.bitpix = 32; precision = 'int32';
case 16,
hdr.dime.bitpix = 32; precision = 'float32';
case 32,
hdr.dime.bitpix = 64; precision = 'float32';
case 64,
hdr.dime.bitpix = 64; precision = 'float64';
case 128,
hdr.dime.bitpix = 24; precision = 'uint8';
case 256
hdr.dime.bitpix = 8; precision = 'int8';
case 511
hdr.dime.bitpix = 96; precision = 'float32';
case 512
hdr.dime.bitpix = 16; precision = 'uint16';
case 768
hdr.dime.bitpix = 32; precision = 'uint32';
case 1024
hdr.dime.bitpix = 64; precision = 'int64';
case 1280
hdr.dime.bitpix = 64; precision = 'uint64';
case 1792,
hdr.dime.bitpix = 128; precision = 'float64';
otherwise
error('This datatype is not supported');
end
hdr.dime.dim(find(hdr.dime.dim < 1)) = 1;
% move pointer to the start of image block
%
switch filetype
case {0, 1}
fseek(fid, 0, 'bof');
case 2
fseek(fid, hdr.dime.vox_offset, 'bof');
end
% Load whole image block for old Analyze format or binary image;
% otherwise, load images that are specified in img_idx, dim5_idx,
% dim6_idx, and dim7_idx
%
% For binary image, we have to read all because pos can not be
% seeked in bit and can not be calculated the way below.
%
if hdr.dime.datatype == 1 | isequal(hdr.dime.dim(5:8),ones(1,4)) | ...
(isempty(img_idx) & isempty(dim5_idx) & isempty(dim6_idx) & isempty(dim7_idx))
% For each frame, precision of value will be read
% in img_siz times, where img_siz is only the
% dimension size of an image, not the byte storage
% size of an image.
%
img_siz = prod(hdr.dime.dim(2:8));
% For complex float32 or complex float64, voxel values
% include [real, imag]
%
if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
img_siz = img_siz * 2;
end
%MPH: For RGB24, voxel values include 3 separate color planes
%
if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
img_siz = img_siz * 3;
end
img = fread(fid, img_siz, sprintf('*%s',precision));
d1 = hdr.dime.dim(2);
d2 = hdr.dime.dim(3);
d3 = hdr.dime.dim(4);
d4 = hdr.dime.dim(5);
d5 = hdr.dime.dim(6);
d6 = hdr.dime.dim(7);
d7 = hdr.dime.dim(8);
if isempty(img_idx)
img_idx = 1:d4;
end
if isempty(dim5_idx)
dim5_idx = 1:d5;
end
if isempty(dim6_idx)
dim6_idx = 1:d6;
end
if isempty(dim7_idx)
dim7_idx = 1:d7;
end
else
d1 = hdr.dime.dim(2);
d2 = hdr.dime.dim(3);
d3 = hdr.dime.dim(4);
d4 = hdr.dime.dim(5);
d5 = hdr.dime.dim(6);
d6 = hdr.dime.dim(7);
d7 = hdr.dime.dim(8);
if isempty(img_idx)
img_idx = 1:d4;
end
if isempty(dim5_idx)
dim5_idx = 1:d5;
end
if isempty(dim6_idx)
dim6_idx = 1:d6;
end
if isempty(dim7_idx)
dim7_idx = 1:d7;
end
% compute size of one image
%
img_siz = prod(hdr.dime.dim(2:4));
% For complex float32 or complex float64, voxel values
% include [real, imag]
%
if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
img_siz = img_siz * 2;
end
%MPH: For RGB24, voxel values include 3 separate color planes
%
if hdr.dime.datatype == 128 | hdr.dime.datatype == 511
img_siz = img_siz * 3;
end
% preallocate img
img = zeros(img_siz, length(img_idx)*length(dim5_idx)*length(dim6_idx)*length(dim7_idx) );
currentIndex = 1;
for i7=1:length(dim7_idx)
for i6=1:length(dim6_idx)
for i5=1:length(dim5_idx)
for t=1:length(img_idx)
% Position is seeked in bytes. To convert dimension size
% to byte storage size, hdr.dime.bitpix/8 will be
% applied.
%
pos = sub2ind([d1 d2 d3 d4 d5 d6 d7], 1, 1, 1, ...
img_idx(t), dim5_idx(i5),dim6_idx(i6),dim7_idx(i7)) -1;
pos = pos * hdr.dime.bitpix/8;
if filetype == 2
fseek(fid, pos + hdr.dime.vox_offset, 'bof');
else
fseek(fid, pos, 'bof');
end
% For each frame, fread will read precision of value
% in img_siz times
%
img(:,currentIndex) = fread(fid, img_siz, sprintf('*%s',precision));
currentIndex = currentIndex +1;
end
end
end
end
end
% For complex float32 or complex float64, voxel values
% include [real, imag]
%
if hdr.dime.datatype == 32 | hdr.dime.datatype == 1792
img = reshape(img, [2, length(img)/2]);
img = complex(img(1,:)', img(2,:)');
end
fclose(fid);
% Update the global min and max values
%
hdr.dime.glmax = double(max(img(:)));
hdr.dime.glmin = double(min(img(:)));
% old_RGB treat RGB slice by slice, now it is treated voxel by voxel
%
if old_RGB & hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
% remove squeeze
img = (reshape(img, [hdr.dime.dim(2:3) 3 hdr.dime.dim(4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
img = permute(img, [1 2 4 3 5 6 7 8]);
elseif hdr.dime.datatype == 128 & hdr.dime.bitpix == 24
% remove squeeze
img = (reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
img = permute(img, [2 3 4 1 5 6 7 8]);
elseif hdr.dime.datatype == 511 & hdr.dime.bitpix == 96
img = double(img(:));
img = single((img - min(img))/(max(img) - min(img)));
% remove squeeze
img = (reshape(img, [3 hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
img = permute(img, [2 3 4 1 5 6 7 8]);
else
% remove squeeze
img = (reshape(img, [hdr.dime.dim(2:4) length(img_idx) length(dim5_idx) length(dim6_idx) length(dim7_idx)]));
end
if ~isempty(img_idx)
hdr.dime.dim(5) = length(img_idx);
end
if ~isempty(dim5_idx)
hdr.dime.dim(6) = length(dim5_idx);
end
if ~isempty(dim6_idx)
hdr.dime.dim(7) = length(dim6_idx);
end
if ~isempty(dim7_idx)
hdr.dime.dim(8) = length(dim7_idx);
end
return % read_image
NIfTI_20140122/load_untouch0_nii_hdr.m 0 → 100644
View file @ ce12fa61
% internal function
% - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
function hdr = load_nii_hdr(fileprefix, machine)
fn = sprintf('%s.hdr',fileprefix);
fid = fopen(fn,'r',machine);
if fid < 0,
msg = sprintf('Cannot open file %s.',fn);
error(msg);
else
fseek(fid,0,'bof');
hdr = read_header(fid);
fclose(fid);
end
return % load_nii_hdr
%---------------------------------------------------------------------
function [ dsr ] = read_header(fid)
% Original header structures
% struct dsr
% {
% struct header_key hk; /* 0 + 40 */
% struct image_dimension dime; /* 40 + 108 */
% struct data_history hist; /* 148 + 200 */
% }; /* total= 348 bytes*/
dsr.hk = header_key(fid);
dsr.dime = image_dimension(fid);
dsr.hist = data_history(fid);
return % read_header
%---------------------------------------------------------------------
function [ hk ] = header_key(fid)
fseek(fid,0,'bof');
% Original header structures
% struct header_key /* header key */
% { /* off + size */
% int sizeof_hdr /* 0 + 4 */
% char data_type[10]; /* 4 + 10 */
% char db_name[18]; /* 14 + 18 */
% int extents; /* 32 + 4 */
% short int session_error; /* 36 + 2 */
% char regular; /* 38 + 1 */
% char hkey_un0; /* 39 + 1 */
% }; /* total=40 bytes */
%
% int sizeof_header Should be 348.
% char regular Must be 'r' to indicate that all images and
% volumes are the same size.
v6 = version;
if str2num(v6(1))<6
directchar = '*char';
else
directchar = 'uchar=>char';
end
hk.sizeof_hdr = fread(fid, 1,'int32')'; % should be 348!
hk.data_type = deblank(fread(fid,10,directchar)');
hk.db_name = deblank(fread(fid,18,directchar)');
hk.extents = fread(fid, 1,'int32')';
hk.session_error = fread(fid, 1,'int16')';
hk.regular = fread(fid, 1,directchar)';
hk.hkey_un0 = fread(fid, 1,directchar)';
return % header_key
%---------------------------------------------------------------------
function [ dime ] = image_dimension(fid)
%struct image_dimension
% { /* off + size */
% short int dim[8]; /* 0 + 16 */
% /*
% dim[0] Number of dimensions in database; usually 4.
% dim[1] Image X dimension; number of *pixels* in an image row.
% dim[2] Image Y dimension; number of *pixel rows* in slice.
% dim[3] Volume Z dimension; number of *slices* in a volume.
% dim[4] Time points; number of volumes in database
% */
% char vox_units[4]; /* 16 + 4 */
% char cal_units[8]; /* 20 + 8 */
% short int unused1; /* 28 + 2 */
% short int datatype; /* 30 + 2 */
% short int bitpix; /* 32 + 2 */
% short int dim_un0; /* 34 + 2 */
% float pixdim[8]; /* 36 + 32 */
% /*
% pixdim[] specifies the voxel dimensions:
% pixdim[1] - voxel width, mm
% pixdim[2] - voxel height, mm
% pixdim[3] - slice thickness, mm
% pixdim[4] - volume timing, in msec
% ..etc
% */
% float vox_offset; /* 68 + 4 */
% float roi_scale; /* 72 + 4 */
% float funused1; /* 76 + 4 */
% float funused2; /* 80 + 4 */
% float cal_max; /* 84 + 4 */
% float cal_min; /* 88 + 4 */
% int compressed; /* 92 + 4 */
% int verified; /* 96 + 4 */
% int glmax; /* 100 + 4 */
% int glmin; /* 104 + 4 */
% }; /* total=108 bytes */
v6 = version;
if str2num(v6(1))<6
directchar = '*char';
else
directchar = 'uchar=>char';
end
dime.dim = fread(fid,8,'int16')';
dime.vox_units = deblank(fread(fid,4,directchar)');
dime.cal_units = deblank(fread(fid,8,directchar)');
dime.unused1 = fread(fid,1,'int16')';
dime.datatype = fread(fid,1,'int16')';
dime.bitpix = fread(fid,1,'int16')';
dime.dim_un0 = fread(fid,1,'int16')';
dime.pixdim = fread(fid,8,'float32')';
dime.vox_offset = fread(fid,1,'float32')';
dime.roi_scale = fread(fid,1,'float32')';
dime.funused1 = fread(fid,1,'float32')';
dime.funused2 = fread(fid,1,'float32')';
dime.cal_max = fread(fid,1,'float32')';
dime.cal_min = fread(fid,1,'float32')';
dime.compressed = fread(fid,1,'int32')';
dime.verified = fread(fid,1,'int32')';
dime.glmax = fread(fid,1,'int32')';
dime.glmin = fread(fid,1,'int32')';
return % image_dimension
%---------------------------------------------------------------------
function [ hist ] = data_history(fid)
%struct data_history
% { /* off + size */
% char descrip[80]; /* 0 + 80 */
% char aux_file[24]; /* 80 + 24 */
% char orient; /* 104 + 1 */
% char originator[10]; /* 105 + 10 */
% char generated[10]; /* 115 + 10 */
% char scannum[10]; /* 125 + 10 */
% char patient_id[10]; /* 135 + 10 */
% char exp_date[10]; /* 145 + 10 */
% char exp_time[10]; /* 155 + 10 */
% char hist_un0[3]; /* 165 + 3 */
% int views /* 168 + 4 */
% int vols_added; /* 172 + 4 */
% int start_field; /* 176 + 4 */
% int field_skip; /* 180 + 4 */
% int omax; /* 184 + 4 */
% int omin; /* 188 + 4 */
% int smax; /* 192 + 4 */
% int smin; /* 196 + 4 */
% }; /* total=200 bytes */
v6 = version;
if str2num(v6(1))<6
directchar = '*char';
else
directchar = 'uchar=>char';
end
hist.descrip = deblank(fread(fid,80,directchar)');
hist.aux_file = deblank(fread(fid,24,directchar)');
hist.orient = fread(fid, 1,'char')';
hist.originator = fread(fid, 5,'int16')';
hist.generated = deblank(fread(fid,10,directchar)');
hist.scannum = deblank(fread(fid,10,directchar)');
hist.patient_id = deblank(fread(fid,10,directchar)');
hist.exp_date = deblank(fread(fid,10,directchar)');
hist.exp_time = deblank(fread(fid,10,directchar)');
hist.hist_un0 = deblank(fread(fid, 3,directchar)');
hist.views = fread(fid, 1,'int32')';
hist.vols_added = fread(fid, 1,'int32')';
hist.start_field = fread(fid, 1,'int32')';
hist.field_skip = fread(fid, 1,'int32')';
hist.omax = fread(fid, 1,'int32')';
hist.omin = fread(fid, 1,'int32')';
hist.smax = fread(fid, 1,'int32')';
hist.smin = fread(fid, 1,'int32')';
return % data_history
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment