1%% Copyright (C) 2012 Adam H Aitkenhead <adamaitkenhead@hotmail.com> 2%% 3%% This program is free software; you can redistribute it and/or modify 4%% it under the terms of the GNU General Public License as published by 5%% the Free Software Foundation; either version 3 of the License, or 6%% (at your option) any later version. 7%% 8%% This program is distributed in the hope that it will be useful, 9%% but WITHOUT ANY WARRANTY; without even the implied warranty of 10%% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11%% GNU General Public License for more details. 12%% 13%% You should have received a copy of the GNU General Public License 14%% along with this program; If not, see <http://www.gnu.org/licenses/>. 15 16%% -*- texinfo -*- 17%% @deftypefn {Function File} analyze75write (@var{filename}, @var{data}, @var{header}) 18%% @deftypefnx {Function File} analyze75write (@var{filename}, @var{data}, @var{x}, @var{y}, @var{z}) 19%% @deftypefnx {Function File} analyze75write (@var{filename}, @var{data}, @var{header}, @var{x}, @var{y}, @var{z}) 20%% Write image data to an Analyze 7.5 file. 21%% 22%% @var{filename} is the path to write the Analyze 7.5 file; @var{data} is 23%% the 3D image data; @var{header} is a structure containing the file 24%% information; @var{x}, @var{y}, @var{z} are lists of the x,y,z coordinates 25%% (in cm) of the data grid. 26%% 27%% @seealso{analyze75info, analyze75read} 28%% @end deftypefn 29 30%% Author: Adam H. Aitkenhead <adamaitkenhead@hotmail.com> 31 32function analyze75write (varargin); 33 34 if (nargin ~= 3 && nargin ~= 5 && nargin ~= 6) 35 print_usage; 36 else 37 filename = varargin{1}; 38 data = varargin{2}; 39 if (nargin ==3) 40 header = varargin{3}; 41 x = ( (0:header.Dimensions(2)-1) - (header.Dimensions(2)-1)/2 ) * header.PixelDimensions(2); 42 y = ( (0:header.Dimensions(1)-1) - (header.Dimensions(1)-1)/2 ) * header.PixelDimensions(1); 43 z = ( (0:header.Dimensions(3)-1) - (header.Dimensions(3)-1)/2 ) * header.PixelDimensions(3); 44 % Convert mm to cm 45 if strncmpi(header.VoxelUnits,'mm',2)==1 46 x = x / 10; 47 y = y / 10; 48 z = z / 10; 49 end 50 elseif (nargin ==5) 51 header = struct; 52 x = varargin{3}; 53 y = varargin{4}; 54 z = varargin{5}; 55 elseif (nargin ==6) 56 header = varargin{3}; 57 x = varargin{4}; 58 y = varargin{5}; 59 z = varargin{6}; 60 end 61 end 62 63 if (~ischar (filename)) 64 error ('analyze75write: `filename'' must be a string.'); 65 end 66 67 %% Strip the filename of the extension 68 fileextH = strfind (filename, '.hdr'); 69 fileextI = strfind (filename, '.img'); 70 if (~isempty (fileextH)) 71 fileprefix = filename(1:fileextH(end)-1); 72 elseif (~isempty (fileextI)) 73 fileprefix = filename(1:fileextI(end)-1); 74 else 75 fileprefix = filename; 76 end 77 78 % Check the byteorder 79 if (nargin == 6) 80 if (isfield (header, 'ByteOrder')) && (any (strcmpi (header.ByteOrder, {'ieee-be', 'b'}))) 81 warning ('analyze75write: No support for big-endian. Please consider submitting a patch. Attempting to write as little-endian'); 82 end 83 end 84 header.ByteOrder = 'ieee-le'; 85 86 %% Rearrange the data 87 data = permute(data,[2,1,3]); 88 89 % Force uniform slice spacing 90 if (max(unique(diff(z))) - min(unique(diff(z))) > 1e-4) 91 warning ('analyze75write: Data slices must be equally spaced. Attempting to interpolate the data onto equally spaced slices.') 92 [data,z] = interpslices(data,x,y,z); 93 end 94 95 % Check certain header fields 96 header.PixelDimensions = single(10 * [x(2)-x(1),y(2)-y(1),z(2)-z(1)]); 97 header.Dimensions = int16([size(data),1]); 98 header.HdrFileSize = int32(348); 99 header.VoxelUnits = 'mm '; 100 header.GlobalMin = int32(min(data(:))); 101 header.GlobalMax = int32(max(data(:))); 102 header.FileName = [fileprefix,'.hdr']; 103 104 % Descriptor: Generate a string containing the coordinates of the first voxel (stored in the header for information only) 105 origintext = ['Coordinates of first voxel: ',num2str( 10* [x(1),y(1),z(1)] ,' %07.2f' )]; 106 header.Descriptor = sprintf('%-80s',origintext); 107 108 % Determine the type of data 109 if (isa(data,'int16')) 110 header.ImgDataType = 'DT_SIGNED_SHORT'; 111 header.BitDepth = int16(16); 112 DataTypeLabel = int16(4); 113 DataTypeString = 'int16'; 114 elseif (isa(data,'int32')) 115 header.ImgDataType = 'DT_SIGNED_INT'; 116 header.BitDepth = int16(32); 117 DataTypeLabel = int16(8); 118 DataTypeString = 'int32'; 119 elseif (isa(data,'single')) 120 header.ImgDataType = 'DT_FLOAT'; 121 header.BitDepth = int16(32); 122 DataTypeLabel = int16(16); 123 DataTypeString = 'single'; 124 elseif (isa(data,'double')) 125 header.ImgDataType = 'DT_DOUBLE'; 126 header.BitDepth = int16(64); 127 DataTypeLabel = int16(64); 128 DataTypeString = 'double'; 129 end 130 131 % Write the .hdr file 132 fidH = fopen([fileprefix,'.hdr'],'w'); 133 134 fwrite(fidH,header.HdrFileSize,'int32',header.ByteOrder); % HdrFileSize 135 fwrite(fidH,repmat(' ',1,10),'char',header.ByteOrder); % HdrDataType 136 fwrite(fidH,repmat(' ',1,18),'char',header.ByteOrder); % DatabaseName 137 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % Extents 138 fwrite(fidH,zeros(1,1,'int16'),'int16',header.ByteOrder); % SessionError 139 fwrite(fidH,'r','char',header.ByteOrder); % Regular 140 fwrite(fidH,repmat(' ',1,1),'char',header.ByteOrder); % unused 141 fwrite(fidH,zeros(1,1,'int16'),'int16',header.ByteOrder); % unused 142 fwrite(fidH,header.Dimensions(1),'int16',header.ByteOrder); % Dimensions(1) 143 fwrite(fidH,header.Dimensions(2),'int16',header.ByteOrder); % Dimensions(2) 144 fwrite(fidH,header.Dimensions(3),'int16',header.ByteOrder); % Dimensions(3) 145 fwrite(fidH,header.Dimensions(4),'int16',header.ByteOrder); % Dimensions(4) 146 fwrite(fidH,zeros(1,3,'int16'),'int16',header.ByteOrder); % unused 147 fwrite(fidH,header.VoxelUnits,'char',header.ByteOrder); % VoxelUnits 148 149 if (isfield(header,'CalibrationUnits')) % CalibrationUnits 150 fwrite(fidH,sprintf('%-8s',header.CalibrationUnits(1:min([8,numel(header.CalibrationUnits)]))),'char',header.ByteOrder); 151 else 152 fwrite(fidH,repmat(' ',1,8),'char',header.ByteOrder); 153 end 154 155 fwrite(fidH,zeros(1,1,'int16'),'int16',header.ByteOrder); % unused 156 fwrite(fidH,DataTypeLabel,'int16',header.ByteOrder); % ImgDataType 157 fwrite(fidH,header.BitDepth,'int16',header.ByteOrder); % BitDepth 158 fwrite(fidH,zeros(1,1,'int16'),'int16',header.ByteOrder); % unused 159 fwrite(fidH,zeros(1,1,'single'),'float',header.ByteOrder); % unused 160 fwrite(fidH,header.PixelDimensions(1),'float',header.ByteOrder); % PixelDimensions(1) 161 fwrite(fidH,header.PixelDimensions(2),'float',header.ByteOrder); % PixelDimensions(2) 162 fwrite(fidH,header.PixelDimensions(3),'float',header.ByteOrder); % PixelDimensions(3) 163 fwrite(fidH,zeros(1,4,'single'),'float',header.ByteOrder); % unused 164 fwrite(fidH,zeros(1,1,'single'),'float',header.ByteOrder); % VoxelOffset 165 fwrite(fidH,zeros(1,3,'single'),'float',header.ByteOrder); % unused 166 fwrite(fidH,zeros(1,1,'single'),'float',header.ByteOrder); % CalibrationMax 167 fwrite(fidH,zeros(1,1,'single'),'float',header.ByteOrder); % CalibrationMin 168 fwrite(fidH,zeros(1,1,'single'),'float',header.ByteOrder); % Compressed 169 fwrite(fidH,zeros(1,1,'single'),'float',header.ByteOrder); % Verified 170 fwrite(fidH,header.GlobalMax,'int32',header.ByteOrder); % GlobalMax 171 fwrite(fidH,header.GlobalMin,'int32',header.ByteOrder); % GlobalMin 172 fwrite(fidH,header.Descriptor,'char',header.ByteOrder); % Descriptor 173 fwrite(fidH,'none ','char',header.ByteOrder); % AuxFile 174 fwrite(fidH,repmat(' ',1,1),'char',header.ByteOrder); % Orientation 175 fwrite(fidH,repmat(' ',1,10),'char',header.ByteOrder); % Originator 176 fwrite(fidH,repmat(' ',1,10),'char',header.ByteOrder); % Generated 177 fwrite(fidH,repmat(' ',1,10),'char',header.ByteOrder); % Scannumber 178 179 if (isfield(header,'PatientID')) % PatientID 180 fwrite(fidH,sprintf('%-10s',header.PatientID(1:min([10,numel(header.PatientID)]))),'char',header.ByteOrder); 181 else 182 fwrite(fidH,repmat(' ',1,10),'char',header.ByteOrder); 183 end 184 185 if (isfield(header,'ExposureDate')) % ExposureDate 186 fwrite(fidH,sprintf('%-10s',header.ExposureDate(1:min([10,numel(header.ExposureDate)]))),'char',header.ByteOrder); 187 elseif (isfield(header,'StudyDate')) 188 fwrite(fidH,sprintf('%-10s',header.StudyDate(1:min([10,numel(header.StudyDate)]))),'char',header.ByteOrder); 189 else 190 fwrite(fidH,repmat(' ',1,10),'char',header.ByteOrder); 191 end 192 193 if (isfield(header,'ExposureTime')) % ExposureTime 194 fwrite(fidH,sprintf('%-10s',header.ExposureTime(1:min([10,numel(header.ExposureTime)]))),'char',header.ByteOrder); 195 elseif (isfield(header,'StudyTime')) 196 fwrite(fidH,sprintf('%-10s',header.StudyTime(1:min([10,numel(header.StudyTime)]))),'char',header.ByteOrder); 197 else 198 fwrite(fidH,repmat(' ',1,10),'char',header.ByteOrder); 199 end 200 201 fwrite(fidH,repmat(' ',1,3),'char',header.ByteOrder); % unused 202 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % Views 203 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % VolumesAdded 204 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % StartField 205 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % FieldSkip 206 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % OMax 207 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % OMin 208 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % SMax 209 fwrite(fidH,zeros(1,1,'int32'),'int32',header.ByteOrder); % SMin 210 211 fclose(fidH); 212 213 % Write the .img file 214 fidI = fopen([fileprefix,'.img'],'w'); 215 fwrite(fidI,data,DataTypeString,header.ByteOrder); 216 fclose(fidI); 217 218end 219 220 221function [datanew,znew] = interpslices(data,x,y,z) 222 znew = z(1):min(diff(z)):z(end); 223 datanew = zeros(numel(x),numel(y),numel(znew),class(data)); 224 for loopN = 1:numel(znew) 225 datanew(:,:,loopN) = interpn(x,y,z,data,x,y,znew(loopN)); 226 end 227end 228