Code covered by the BSD License  

Highlights from
Tools for NIfTI and ANALYZE image

image thumbnail

Tools for NIfTI and ANALYZE image

by

Jimmy Shen (view profile)

 

23 Oct 2005 (Updated )

Load, save, make, reslice, view (and edit) both NIfTI and ANALYZE data on any platform

save_untouch_slice(slice, filename, slice_idx, img_idx, dim5_idx, dim6_idx, dim7_idx)
%  Save back to the original image with a portion of slices that was
%  loaded by "load_untouch_nii". You can process those slices matrix
%  in any way, as long as their dimension is not altered.
%
%  Usage: save_untouch_slice(slice, filename, ...
%		slice_idx, [img_idx], [dim5_idx], [dim6_idx], [dim7_idx])
%
%  slice  -  a portion of slices that was loaded by "load_untouch_nii".
%	This should be a numeric matrix (i.e. only the .img field in the
%	loaded structure)
%
%  filename  - 	NIfTI or ANALYZE file name.
%
%  slice_idx (depending on slice size)  -  a numerical array of image
%	slice indices, which should be the same as that you entered
%	in "load_untouch_nii" command.
%
%  img_idx (depending on slice size)  -  a numerical array of image
%	volume indices, which should be the same as that you entered
%	in "load_untouch_nii" command.
%
%  dim5_idx (depending on slice size)  -  a numerical array of 5th 
%	dimension indices, which should be the same as that you entered
%	in "load_untouch_nii" command.
%
%  dim6_idx (depending on slice size)  -  a numerical array of 6th 
%	dimension indices, which should be the same as that you entered
%	in "load_untouch_nii" command.
%
%  dim7_idx (depending on slice size)  -  a numerical array of 7th 
%	dimension indices, which should be the same as that you entered
%	in "load_untouch_nii" command.
%
%  Example:
%	nii = load_nii('avg152T1_LR_nifti.nii');
%	save_nii(nii, 'test.nii');
%	view_nii(nii);
%	nii = load_untouch_nii('test.nii','','','','','',[40 51:53]);
%	nii.img = ones(91,109,4)*122;
%	save_untouch_slice(nii.img, 'test.nii', [40 51:52]);
%	nii = load_nii('test.nii');
%	view_nii(nii);
%
%  - Jimmy Shen (jimmy@rotman-baycrest.on.ca)
%
function save_untouch_slice(slice, filename, slice_idx, img_idx, dim5_idx, dim6_idx, dim7_idx)

   if ~exist('slice','var') | ~isnumeric(slice)
      msg = [char(10) '"slice" argument should be a portion of slices that was loaded' char(10)];
      msg = [msg 'by "load_untouch_nii.m". This should be a numeric matrix (i.e.' char(10)];
      msg = [msg 'only the .img field in the loaded structure).'];
      error(msg);
   end

   if ~exist('filename','var') | ~exist(filename,'file')
      error('In order to save back, original NIfTI or ANALYZE file must exist.');
   end

   if ~exist('slice_idx','var') | isempty(slice_idx) | ~isequal(size(slice,3),length(slice_idx))
      msg = [char(10) '"slice_idx" is a numerical array of image slice indices, which' char(10)];
      msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
      msg = [msg 'command.'];
      error(msg);
   end

   if ~exist('img_idx','var') | isempty(img_idx)
      img_idx = [];

      if ~isequal(size(slice,4),1)
         msg = [char(10) '"img_idx" is a numerical array of image volume indices, which' char(10)];
         msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
         msg = [msg 'command.'];
         error(msg);
      end
   elseif ~isequal(size(slice,4),length(img_idx))
      msg = [char(10) '"img_idx" is a numerical array of image volume indices, which' char(10)];
      msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
      msg = [msg 'command.'];
      error(msg);
   end

   if ~exist('dim5_idx','var') | isempty(dim5_idx)
      dim5_idx = [];

      if ~isequal(size(slice,5),1)
         msg = [char(10) '"dim5_idx" is a numerical array of 5th dimension indices, which' char(10)];
         msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
         msg = [msg 'command.'];
         error(msg);
      end
   elseif ~isequal(size(slice,5),length(img_idx))
      msg = [char(10) '"img_idx" is a numerical array of 5th dimension indices, which' char(10)];
      msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
      msg = [msg 'command.'];
      error(msg);
   end

   if ~exist('dim6_idx','var') | isempty(dim6_idx)
      dim6_idx = [];

      if ~isequal(size(slice,6),1)
         msg = [char(10) '"dim6_idx" is a numerical array of 6th dimension indices, which' char(10)];
         msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
         msg = [msg 'command.'];
         error(msg);
      end
   elseif ~isequal(size(slice,6),length(img_idx))
      msg = [char(10) '"img_idx" is a numerical array of 6th dimension indices, which' char(10)];
      msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
      msg = [msg 'command.'];
      error(msg);
   end

   if ~exist('dim7_idx','var') | isempty(dim7_idx)
      dim7_idx = [];

      if ~isequal(size(slice,7),1)
         msg = [char(10) '"dim7_idx" is a numerical array of 7th dimension indices, which' char(10)];
         msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
         msg = [msg 'command.'];
         error(msg);
      end
   elseif ~isequal(size(slice,7),length(img_idx))
      msg = [char(10) '"img_idx" is a numerical array of 7th dimension indices, which' char(10)];
      msg = [msg 'should be the same as that you entered in "load_untouch_nii.m"' char(10)];
      msg = [msg 'command.'];
      error(msg);
   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);

   if nii.filetype == 0
      nii.hdr = load_untouch0_nii_hdr(nii.fileprefix,nii.machine);
   else
      nii.hdr = load_untouch_nii_hdr(nii.fileprefix,nii.machine,nii.filetype);
   end


   %  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

   [p,f] = fileparts(filename);
   fileprefix = fullfile(p, f);
%   fileprefix = nii.fileprefix;
   filetype = nii.filetype;

   if ~isequal( nii.hdr.dime.dim(2:3), [size(slice,1),size(slice,2)] )
      msg = [char(10) 'The first two dimensions of slice matrix should be the same as' char(10)];
      msg = [msg 'the first two dimensions of image loaded by "load_untouch_nii".'];
      error(msg);
   end


   %  Save the dataset body
   %
   save_untouch_slice_img(slice, nii.hdr, filetype, fileprefix, ...
	nii.machine, slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx);

   %  gzip output file if requested
   %
   if exist('gzFileName', 'var')
      [p,f] = fileparts(gzFileName);

      if filetype == 1
         gzip([fileprefix, '.img']);
         delete([fileprefix, '.img']);
         movefile([fileprefix, '.img.gz']);
         gzip([fileprefix, '.hdr']);
         delete([fileprefix, '.hdr']);
         movefile([fileprefix, '.hdr.gz']);
      elseif filetype == 2
         gzip([fileprefix, '.nii']);
         delete([fileprefix, '.nii']);
         movefile([fileprefix, '.nii.gz']);
      end;

      rmdir(tmpDir,'s');
   end;

   return					% save_untouch_slice


%--------------------------------------------------------------------------
function save_untouch_slice_img(slice,hdr,filetype,fileprefix,machine,slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx)

   if ~exist('hdr','var') | ~exist('filetype','var') | ~exist('fileprefix','var') | ~exist('machine','var')
      error('Usage: save_untouch_slice_img(slice,hdr,filetype,fileprefix,machine,slice_idx,[img_idx],[dim5_idx],[dim6_idx],[dim7_idx]);');
   end

   if ~exist('slice_idx','var') | isempty(slice_idx) | hdr.dime.dim(4)<1
      slice_idx = [];
   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

   %  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

   %  check slice_idx
   %
   if ~isempty(slice_idx) & ~isnumeric(slice_idx)
      error('"slice_idx" should be a numerical array.');
   end

   if length(unique(slice_idx)) ~= length(slice_idx)
      error('Duplicate index in "slice_idx"');
   end

   if ~isempty(slice_idx) & (min(slice_idx) < 1 | max(slice_idx) > hdr.dime.dim(4))
      max_range = hdr.dime.dim(4);

      if max_range == 1
         error(['"slice_idx" should be 1.']);
      else
         range = ['1 ' num2str(max_range)];
         error(['"slice_idx" should be an integer within the range of [' range '].']);
      end
   end

   write_image(slice,hdr,filetype,fileprefix,machine,slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx);

   return					% save_untouch_slice_img


%---------------------------------------------------------------------
function write_image(slice,hdr,filetype,fileprefix,machine,slice_idx,img_idx,dim5_idx,dim6_idx,dim7_idx)

   if filetype == 2
      fid = fopen(sprintf('%s.nii',fileprefix),'r+');

      if fid < 0,
         msg = sprintf('Cannot open file %s.nii.',fileprefix);
         error(msg);
      end
   else
      fid = fopen(sprintf('%s.img',fileprefix),'r+');

      if fid < 0,
         msg = sprintf('Cannot open file %s.img.',fileprefix);
         error(msg);
      end
   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   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  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';
   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

   if hdr.dime.datatype == 1 | isequal(hdr.dime.dim(4:8),ones(1,5)) | ...
	(isempty(img_idx) & isempty(dim5_idx) & isempty(dim6_idx) & isempty(dim7_idx) & isempty(slice_idx))

      msg = [char(10) char(10) '   "save_untouch_slice" is used to save back to the original image a' char(10)];
      msg = [msg '   portion of slices that were loaded by "load_untouch_nii". You can' char(10)];
      msg = [msg '   process those slices matrix in any way, as long as their dimension' char(10)];
      msg = [msg '   is not changed.'];
      error(msg);
   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(slice_idx)
         slice_idx = 1:d3;
      end

      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
      
      %ROMAN: begin
      roman = 1;
      if(roman)

         %  compute size of one slice
         %
         img_siz = prod(hdr.dime.dim(2:3));

         %  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

      end; %if(roman)
      % ROMAN: end

      for i7=1:length(dim7_idx)
         for i6=1:length(dim6_idx)
            for i5=1:length(dim5_idx)
               for t=1:length(img_idx)
               for s=1:length(slice_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, slice_idx(s), ...
			                    img_idx(t), dim5_idx(i5),dim6_idx(i6),dim7_idx(i7)) -1;
                  pos = pos * hdr.dime.bitpix/8;

                  % ROMAN: begin
                  if(roman)
                      % do nothing
                  else
                     img_siz = prod(hdr.dime.dim(2:3));

                     %  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
                  end; % if (roman)
                  % ROMAN: end
         
                  if filetype == 2
                     fseek(fid, pos + hdr.dime.vox_offset, 'bof');
                  else
                     fseek(fid, pos, 'bof');
                  end

                  %  For each frame, fwrite will write precision of value
                  %  in img_siz times
                  %
                  fwrite(fid, slice(:,:,s,t,i5,i6,i7), sprintf('*%s',precision));
                  
               end
               end
            end
         end
      end
   end

   fclose(fid);

   return						% write_image

Contact us