matRad_writeNRRD

 matRad NRRD writer
 
 call
   matRad_writeNRRD(filename,cube,datatype,...
                    additionalFields,additionalKeyValuePairs)

 input
   filename:   full output path, including the nrrd extension
   cube:       cube that is to be written
   metadata:   struct of metadata. Writer will wrap the existing metadata 
               to nrrd standard-specific fields [1].

 References
   [1] http://teem.sourceforge.net/nrrd/format.html5

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 Copyright 2015 the matRad development team. 
 
 This file is part of the matRad project. It is subject to the license 
 terms in the LICENSE file found in the top-level directory of this 
 distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part 
 of the matRad project, including this file, may be copied, modified, 
 propagated, or distributed except according to the terms contained in the 
 LICENSE file.

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

0001 function matRad_writeNRRD(filename,cube,metadata)
0002 % matRad NRRD writer
0003 %
0004 % call
0005 %   matRad_writeNRRD(filename,cube,datatype,...
0006 %                    additionalFields,additionalKeyValuePairs)
0007 %
0008 % input
0009 %   filename:   full output path, including the nrrd extension
0010 %   cube:       cube that is to be written
0011 %   metadata:   struct of metadata. Writer will wrap the existing metadata
0012 %               to nrrd standard-specific fields [1].
0013 %
0014 % References
0015 %   [1] http://teem.sourceforge.net/nrrd/format.html5
0016 %
0017 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0018 %
0019 % Copyright 2015 the matRad development team.
0020 %
0021 % This file is part of the matRad project. It is subject to the license
0022 % terms in the LICENSE file found in the top-level directory of this
0023 % distribution and at https://github.com/e0404/matRad/LICENSES.txt. No part
0024 % of the matRad project, including this file, may be copied, modified,
0025 % propagated, or distributed except according to the terms contained in the
0026 % LICENSE file.
0027 %
0028 % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0029     
0030 %% Setup Header
0031 %The NRRD header description line
0032 version = 5;
0033 nrrdVersionStringPrefix = 'NRRD000';
0034 nrrdVersionString = [nrrdVersionStringPrefix num2str(version)];
0035 
0036 header = sprintf('%s\n',nrrdVersionString);
0037 
0038 %add matRad specific comment
0039 header = header_addComment(header,'Created With matRad - An open source multi-modality radiation treatment planning sytem');
0040 
0041 %Add Datatype field
0042 header = header_addField(header,'type',matlabTypeToNRRD(metadata.datatype));
0043 %NRRD wants to know how the bytes are aligned for datatypes > 1 byte
0044 %(little endian vs big endian)
0045 try
0046     z = zeros(1,metadata.datatype);
0047     varInfo = whos('z');
0048     if varInfo.bytes > 1
0049         [~,~,endian] = computer;
0050         switch endian
0051             case 'L'
0052                 header = header_addField(header,'endian','little');
0053             case 'B'
0054                 header = header_addField(header,'endian','big');
0055             otherwise
0056                 error('Unknown endian!');
0057         end
0058     end
0059 catch
0060     error(['Unknown datatype: ' metadata.datatype']);
0061 end
0062 
0063 
0064 %Check if compression was defined
0065 if ~isfield(metadata,'compress')
0066     metadata.compress = false;
0067 end  
0068 %If we use compression, try to compress, otherwise use raw data
0069 if metadata.compress
0070     try
0071         %We write a temporary file to gzip with matlab
0072         %It might be also possible to do this with java directly as stream
0073         tmpRawFile = [tempname() '.bin'];
0074         fTmp = fopen(tmpRawFile, 'wb');
0075         
0076         fwrite(fTmp,cube,metadata.datatype);
0077         fclose(fTmp);
0078         
0079         %Now zip the file
0080         fNameZip = gzip(tmpRawFile);
0081         fNameZip = fNameZip{1};
0082         
0083         %We need this to determine the file size
0084         fDir = dir(fNameZip);
0085         
0086         %Read the zipped data
0087         fTmp = fopen(fNameZip, 'rb');
0088         fileData = fread(fTmp,fDir.bytes,'uint8');
0089         fclose(fTmp);
0090         %Delete temporary files
0091         delete(tmpRawFile,fNameZip);
0092         
0093         %Set the datatype to binary for writing the zipped data
0094         metadata.datatype = 'uint8';
0095         header = header_addField(header,'encoding','gzip');
0096     catch
0097         warn('Could not open temporary file, writing without compression!');
0098         header = header_addField(header,'encoding','raw');
0099         fileData = cube;
0100     end
0101 else
0102     header = header_addField(header,'encoding','raw');
0103     fileData = cube;
0104 end
0105 
0106 
0107 %Dimensionality
0108 cubeDim = size(cube);
0109 header = header_addField(header,'dimension',num2str(numel(cubeDim)));
0110 strCubeDim = vec2str(cubeDim);
0111 header = header_addField(header,'sizes',strCubeDim);
0112 
0113 %Coordinate system - optional
0114 if isfield(metadata,'coordinateSystem')
0115     header = header_addField(header,'space',metadata.coordinateSystem);
0116 end
0117 
0118 %Reference Point - optional
0119 if isfield(metadata,'imageOrigin')
0120     vecOrigin = matVec2nrrdVec(metadata.imageOrigin);
0121     header = header_addField(header,'space origin',vecOrigin);
0122 end
0123 
0124 %Axis permutation
0125 if isfield(metadata,'resolution')
0126     %If we have more space information, write the space axis
0127     if isfield(metadata,'axisPermutation')
0128         strDirection = '';
0129         for dim = 1:numel(cubeDim)
0130             dimVec = zeros(1,numel(cubeDim));
0131             dimVec(dim) = 1;
0132             dimVec = dimVec(metadata.axisPermutation);
0133             dimVec = dimVec .* metadata.resolution;
0134             strDirection = [strDirection matVec2nrrdVec(dimVec) ' '];
0135         end
0136         header = header_addField(header,'space directions',strDirection);
0137     %Otherwise only write resolution
0138     else
0139         strSpacings = vec2str(metadata.resolution);
0140         header = header_addField(header,'spacings',strSpacings);
0141     end
0142 end       
0143 
0144 %Additional key-value-pairs if present
0145 %Will be added in a later release
0146 %{
0147 if nargin > 3
0148     keys = fieldnames(additionalKeyValuePairs);
0149     for f = 1:numel(fields)
0150         key = keys{f};
0151         value = additionalKeyValuePairs.(keys{f});
0152         
0153         %Format the description to fit
0154         if isvector(description)
0155             value = mat2str(description);
0156             value = description(2:end-1);
0157         elseif isnumeric(description)
0158             value = num2str(description);
0159         else 
0160      
0161         end
0162         
0163         if and(~ischar(value),~iscellstr(value))
0164             warning(['Cannot convert value for key ' key ' to string, key-value-pair will be ignored!']);
0165             continue;
0166         end
0167         
0168         header = header_addKeyValuePair(header,key,value);
0169     end;
0170 end
0171 %}
0172 
0173 %% Write File
0174 try
0175     %Write Header to file with the separating blank line
0176     fileHandle = fopen(filename,'w');
0177     fprintf(fileHandle,'%s\n',header);
0178     
0179     %Write raw or compressed data to file
0180     fwrite(fileHandle,fileData,metadata.datatype);
0181     fclose(fileHandle);
0182     
0183 catch MExc
0184     fclose('all');
0185     error(sprintf('File %s could not be written!\n%s',filename,getReport(MExc)));
0186 end
0187 
0188 end
0189 
0190 %Used to add comments to the header
0191 function newHeader = header_addComment(header,comment)
0192     newHeader = sprintf('%s# %s\n',header,comment);
0193 end
0194 
0195 %Used to add fields to the header
0196 function newHeader = header_addField(header,fieldName,fieldDescription)
0197     newHeader = sprintf('%s%s: %s\n',header,fieldName,fieldDescription);
0198 end
0199 
0200 %Used to add key-value pairs to the header
0201 function newHeader = header_addKeyValuePair(header,key,value)
0202     value = strrep(value,':=','');
0203     newHeader = sprintf('%s%s:=%s\n',header,key,value);
0204 end
0205 
0206 %Used to format dimenional info according to nrrd specification
0207 function strOutput = vec2str(V)
0208     strOutput = mat2str(V);
0209     strOutput = strOutput(2:end-1);
0210 end
0211 
0212 %Used to format a matlab vector to nrrd vector with brackets
0213 function strOutput = matVec2nrrdVec(V)
0214     strOutput = '(';
0215     for v = 1:numel(V)
0216         strOutput = [strOutput num2str(V(v)) ','];
0217     end
0218     strOutput = [strOutput(1:end-1) ')'];
0219 end
0220 
0221 %Used to map matlab definitions of datatypes to the NRRD specs
0222 %seems to be only single has to be converted to float
0223 function newType = matlabTypeToNRRD(datatype)
0224 switch datatype
0225     case 'single'
0226         newType = 'float';
0227     otherwise
0228         newType = datatype;
0229 end
0230 end
0231