spm_file_merge.m 5.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180
  1. function V4 = spm_file_merge(V,fname,dt,RT)
  2. % Concatenate 3D volumes into a single 4D volume
  3. % FUNCTION V4 = spm_file_merge(V,fname,dt)
  4. % V - images to concatenate (char array or spm_vol struct)
  5. % fname - filename for output 4D volume [defaults: '4D.nii']
  6. % Unless explicit, output folder is the one containing first image
  7. % dt - datatype (see spm_type) [defaults: 0]
  8. % 0 means same datatype than first input volume
  9. % RT - Interscan interval {seconds} [defaults: NaN]
  10. %
  11. % V4 - spm_vol struct of the 4D volume
  12. %__________________________________________________________________________
  13. %
  14. % For integer datatypes, the file scale factor is chosen as to maximise
  15. % the range of admissible values. This may lead to quantization error
  16. % differences between the input and output images values.
  17. %__________________________________________________________________________
  18. % Copyright (C) 2009-2018 Wellcome Trust Centre for Neuroimaging
  19. % John Ashburner
  20. % $Id: spm_file_merge.m 7354 2018-06-22 10:44:22Z guillaume $
  21. %-Input: V
  22. %--------------------------------------------------------------------------
  23. if ~nargin
  24. [V,sts] = spm_select([1 Inf],'image','Select images to concatenate');
  25. if ~sts, return; end
  26. end
  27. if ischar(V)
  28. V = spm_vol(V);
  29. elseif iscellstr(V)
  30. V = spm_vol(char(V));
  31. end
  32. ind = cat(1,V.n);
  33. N = cat(1,V.private);
  34. %-Input: fname
  35. %--------------------------------------------------------------------------
  36. cwd = fileparts(V(1).fname);
  37. if isempty(cwd), cwd = pwd; end
  38. if nargin < 2
  39. fname = fullfile(cwd,'4D.nii');
  40. else
  41. [p,n,e] = fileparts(fname);
  42. if isempty(e), fname = [fname '.nii']; end
  43. if isempty(p), fname = fullfile(cwd, fname); end
  44. end
  45. %-Input: dtype
  46. %--------------------------------------------------------------------------
  47. if nargin < 3
  48. dt = 0;
  49. end
  50. if dt == 0
  51. dt = V(1).dt(1);
  52. end
  53. %-Input: RT
  54. %--------------------------------------------------------------------------
  55. if nargin < 4
  56. RT = NaN;
  57. end
  58. if isnan(RT) && ...
  59. isfield(V(1).private,'timing') && isfield(V(1).private.timing,'tspace')
  60. RT = V(1).private.timing.tspace;
  61. end
  62. %-Set scalefactors and offsets
  63. %==========================================================================
  64. d = cat(1,V.dt); d = d(:,1);
  65. s = cat(2,V.pinfo); s = s(1,:);
  66. o = cat(2,V.pinfo); o = o(2,:);
  67. %-Reuse parameters of input images if same scalefactor, offset and datatype
  68. %--------------------------------------------------------------------------
  69. if length(unique(s)) == 1 && length(unique(o)) == 1 ...
  70. && length(unique(d)) == 1 && d(1) == dt
  71. sf = V(1).pinfo(1);
  72. off = V(1).pinfo(2);
  73. else
  74. dmx = spm_type(dt,'maxval');
  75. dmn = spm_type(dt,'minval');
  76. %-Integer datatypes: scale to min/max of overall data
  77. %----------------------------------------------------------------------
  78. if isfinite(dmx)
  79. spm_progress_bar('Init',numel(V),'Computing scale factor','Volumes Complete');
  80. mx = -Inf;
  81. mn = Inf;
  82. for i=1:numel(V)
  83. dat = V(i).private.dat(:,:,:,ind(i,1),ind(i,2));
  84. dat = dat(isfinite(dat));
  85. mx = max(mx,max(dat(:)));
  86. mn = min(mn,min(dat(:)));
  87. spm_progress_bar('Set',i);
  88. end
  89. spm_progress_bar('Clear');
  90. if isempty(mx), mx = 0; end
  91. if isempty(mn), mn = 0; end
  92. if mx~=mn
  93. if dmn < 0
  94. sf = max(mx/dmx,-mn/dmn);
  95. else
  96. sf = mx/dmx;
  97. end
  98. off = 0;
  99. else
  100. sf = mx/dmx;
  101. off = 0;
  102. end
  103. %-floating precison: no scaling
  104. %----------------------------------------------------------------------
  105. else
  106. sf = 1;
  107. off = 0;
  108. end
  109. end
  110. %-Create and write 4D volume image
  111. %==========================================================================
  112. spm_unlink(fname);
  113. %-Create NifTI header
  114. %--------------------------------------------------------------------------
  115. ni = nifti;
  116. ni.dat = file_array(fname,...
  117. [V(1).dim numel(V)],...
  118. [dt spm_platform('bigend')],...
  119. 0,...
  120. sf,...
  121. off);
  122. ni.mat = N(1).mat;
  123. ni.mat0 = N(1).mat;
  124. ni.descrip = '4D image';
  125. if ~isnan(RT)
  126. ni.timing = struct('toffset',0, 'tspace',RT);
  127. end
  128. create(ni);
  129. %-Write 4D data
  130. %--------------------------------------------------------------------------
  131. spm_progress_bar('Init',size(ni.dat,4),'Saving 4D image','Volumes Complete');
  132. for i=1:size(ni.dat,4)
  133. ni.dat(:,:,:,i) = N(i).dat(:,:,:,ind(i,1),ind(i,2));
  134. spm_get_space([ni.dat.fname ',' num2str(i)], V(i).mat);
  135. spm_progress_bar('Set',i);
  136. end
  137. spm_progress_bar('Clear');
  138. %-Fix ?form_code in header (mat_intent is changed by spm_get_space above)
  139. %--------------------------------------------------------------------------
  140. ni = nifti(fname);
  141. ni.mat_intent = N(1).mat_intent;
  142. ni.mat0_intent = N(1).mat0_intent;
  143. create(ni);
  144. %-Remove .mat file if present and not necessary
  145. %--------------------------------------------------------------------------
  146. matfname = spm_file(fname,'ext','mat');
  147. if spm_existfile(matfname)
  148. M = load(matfname);
  149. if isequal(fieldnames(M),{'mat'}) % contains only 'mat'
  150. if sum(sum(M.mat(:,:,1).^2))==0
  151. M.mat(:,:,1) = N(1).mat;
  152. end
  153. if sum(sum(diff(M.mat,1,3).^2))<1e-8
  154. spm_unlink(matfname);
  155. end
  156. end
  157. end
  158. %-Return spm_vol structure
  159. %--------------------------------------------------------------------------
  160. if nargout
  161. V4 = spm_vol(fname);
  162. end