Description of matRad

0001 function [cube, metadata] = matRad_readNRRD(filename)
0002 % matRad NRRD reader
0003 %
0004 % call
0005 %   [cube, metadata] = matRad_readNRRD(filename)
0006 %
0007 % input
0008 %   filename:   full path to nrrd file
0009 %
0010 % output
0011 %   cube:       the read cube
0012 %   metadata:   metadata from header information
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 % Open file.
0031 hFile = fopen(filename, 'r');
0032 if hFile <= 0
0033     error('Could not open NRRD file!');
0034 end
0035 cleaner = onCleanup(@() fclose(hFile));
0036 
0037 %% Determine NRRD Version
0038 nrrdLine = fgetl(hFile);
0039 regTokens = regexp(nrrdLine,'NRRD00(?:\.)?0([1-5])','tokens');
0040 if isempty(regTokens)
0041     error('Invalid Header line! Could not identify NRRD version!');
0042 end
0043 nrrdVersion = str2num(regTokens{1}{1});
0044 
0045 if nrrdVersion > 5
0046     error('NRRD version > 5 not supported!');
0047 end
0048 
0049 %% Read header
0050 nrrdMetaData.fields = cell(0,2);
0051 nrrdMetaData.keys = cell(0,2);
0052 nrrdMetaData.comments = cell(0);
0053 
0054 currentLine = fgetl(hFile);
0055 while ~isempty(currentLine) && ischar(currentLine) %NRRD separates data from header by an empty line
0056     %Check for comment
0057     if isequal(currentLine(1),'#')
0058         nrrdMetaData.comments{end+1} = currentLine;
0059     else
0060         %Parse the line
0061         lineContent = regexp(currentLine, '(.+):(=|\s)(.+)', 'tokens');
0062         if isempty(lineContent)
0063             warning(['Could not parse line: "' lineContent '"']);
0064         elseif isequal(lineContent{1}{2},' ') %space after colon refers to "field"
0065             nrrdMetaData.fields{end+1,1} = lineContent{1}{1}; %Fieldname
0066             nrrdMetaData.fields{end,2} = lineContent{1}{3}; %Information
0067         elseif isequal(lineContent{1}{2},'=') %= after colon refers to key-value pair
0068             nrrdMetaData.keys{end+1,1} = lineContent{1}{1}; %Key
0069             nrrdMetaData.keys{end,2} = lineContent{1}{3}; %Value
0070         else
0071             warning(['Could not parse line: "' lineContent '"']);
0072         end        
0073     end
0074     currentLine = fgetl(hFile);
0075 end
0076 
0077 %% Interpret Headers
0078 %check for the always required data type (and endian if required)
0079 doSwapBytes = false; %If endian differs
0080 typeFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'type'));
0081 if ~isempty(typeFieldIx)
0082     datatype = getMATLABdataType(nrrdMetaData.fields{typeFieldIx,2});
0083     %if size of the datatype is > 1 byte, we need endian information
0084     if ~isequal(datatype(end),'8')
0085         endianFieldIx = find(ismember(nrrdMetaData.fields(:,1),'endian'));
0086         if ~isempty(endianFieldIx)
0087             if ~isequal(nrrdMetaData.fields{endianFieldIx,2},'little') && ~isequal(nrrdMetaData.fields{endianFieldIx,2},'big')
0088                 error(['Datatype is ' datatype ', thus endian information is required but could not be interpreted!']);
0089             end;
0090             %Now we compare the file endian to the system endian
0091             %First acquire system endian
0092             [~,~,endian] = computer();
0093             if isequal(endian,'B')
0094                 endian = 'big';
0095             end;
0096             if isequal(endian,'L')
0097                 endian = 'little';
0098             end;
0099             %now compare to file endian and set flag if appropriate
0100             if ~isequal(endian,nrrdMetaData.fields{endianFieldIx,2})
0101                 doSwapBytes = true;
0102             end
0103         else
0104             error(['Datatype is ' datatype ', thus endian information is required but could not be found!']);
0105         end
0106     end
0107     metadata.datatype = datatype;
0108     
0109 else
0110     error('Could not find required "type" field!');
0111 end
0112 
0113 %Check for the always required image dimension
0114 dimFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'dimension')); 
0115 if ~isempty(dimFieldIx)
0116     [metadata.dimension,success] = str2num(nrrdMetaData.fields{dimFieldIx,2});
0117     if ~success
0118         error('Could not read required dimension field');
0119     end
0120 else
0121     error('Could not find required "dimension" field!');
0122 end
0123 
0124 %Check for size / dim length
0125 sizeFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'sizes')); 
0126 if ~isempty(sizeFieldIx)
0127     sizes = textscan(nrrdMetaData.fields{sizeFieldIx,2},'%d');
0128     if numel(sizes{1}) ~= metadata.dimension || ~all(sizes{1} > 0) 
0129         error('Incorrect size definition!');
0130     end
0131     metadata.cubeDim = sizes{1}';
0132 else
0133     error('Could not find required "dimension" field!');
0134 end
0135 
0136 %Check for resolution
0137 %Here we need either spacings ore space directions
0138 spacingFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'spacings'));
0139 spaceDirFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'space directions'));
0140 if ~isempty(spacingFieldIx)
0141     resolutions = textscan(nrrdMetaData.fields{spacingFieldIx,2},'%f');
0142     if numel(resolutions{1}) ~= metadata.dimension
0143         error('Incorrect spacings definition');
0144     end
0145     metadata.resolution = resolutions{1}';    
0146     
0147     %We have default permutation here (required for mapping to MATLAB
0148     %default order 2 1 3)
0149     metadata.axisPermutation = [1 2 3];
0150     
0151 elseif ~isempty(spaceDirFieldIx)
0152     %space directions are written in vector format
0153     %first create the scanning string
0154     vectorstring = '(';
0155     for c=1:metadata.dimension
0156         vectorstring = [vectorstring '%f,'];
0157     end
0158     vectorstring(end) = ')';
0159     %Get the vectors
0160     vectors = textscan(nrrdMetaData.fields{spaceDirFieldIx,2},vectorstring);
0161     
0162     %At the moment we only support cartesian basis vectors
0163     %this gives us the permutation to align with the MATLAB default
0164     %ordering of 2 1 3
0165     for c=1:metadata.dimension
0166         %check if cartesian basis vector
0167         currentAxis = find(vectors{c});
0168         
0169         if numel(find(vectors{c})) ~= 1
0170             error('Sorry! We currently only support spaces with cartesian basis!'); 
0171         end        
0172         metadata.axisPermutation(c) = currentAxis*sign(vectors{c}(currentAxis));
0173         metadata.resolution(c) = vectors{c}(currentAxis);       
0174     end
0175    
0176 else
0177     warning('No Resolution Information available');
0178 end
0179 
0180 %find the origin if we have one
0181 originFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'space origin'));
0182 if ~isempty(originFieldIx)
0183     %first create the scanning string
0184     vectorstring = '(';
0185     for c=1:metadata.dimension
0186         vectorstring = [vectorstring '%f,'];
0187     end
0188     originVector = textscan(nrrdMetaData.fields{originFieldIx,2},vectorstring);
0189     for c=1:metadata.dimension
0190         metadata.imageOrigin(c) = originVector{c};
0191     end
0192     %metadata.imageOrigin = transpose(metadata.imageOrigin);
0193 end
0194 
0195 %Coordinate system - optional
0196 originFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'space'));
0197 if ~isempty(originFieldIx)
0198     metadata.coordinateSystem = nrrdMetaData.fields{originFieldIx,2};
0199 end
0200     
0201 %Check for separate file
0202 data_fileFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'data file'));
0203 datafileFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'datafile'));
0204 if ~isempty(data_fileFieldIx) || ~isempty(datafileFieldIx)
0205     error('Sorry! We currently do not support detached data files!');
0206     %Proposed workflow:
0207     %check for data file
0208     %close file
0209     %replace file handle with data file
0210     %continue without the read part knowing its a new file
0211 end
0212 
0213 
0214 %% Read Data
0215 %Check for encoding
0216 encodingFieldIx = find(ismember(nrrdMetaData.fields(:,1), 'encoding'));
0217 if isempty(encodingFieldIx)
0218     error('Could not find required "encoding" field!');
0219 end
0220 switch nrrdMetaData.fields{encodingFieldIx,2}
0221     case 'raw'
0222         cube = fread(hFile,prod(metadata.cubeDim), metadata.datatype);
0223     case {'txt','text','ascii'}
0224         cube = cast(fscanf(hFile,'%f'),metadata.datatype);
0225     case 'hex'
0226         error('Sorry: NRRD hex file not yet supported!');
0227     case {'gz','gzip'}
0228         compressedByteArray = fread(hFile, inf, 'uint8');
0229         
0230         javaUnpackSuccessful = false;
0231         
0232         if isempty(javachk('jvm'))
0233             %Unzip with java
0234             try
0235                 javaByteStream = java.io.ByteArrayInputStream(compressedByteArray);
0236                 gzipStream = java.util.zip.GZIPInputStream(javaByteStream);
0237                 javaByteOutputStream = java.io.ByteArrayOutputStream();
0238                 org.apache.commons.io.IOUtils.copy(gzipStream,javaByteOutputStream);
0239                 gzipStream.close();
0240                 cube = typecast(javaByteOutputStream.toByteArray(),metadata.datatype)';
0241                 javaByteStream.close();
0242                 javaByteOutputStream.close();
0243                 javaUnpackSuccessful = true;
0244             catch
0245                 warning('Java unpacking failed... using temporary files!');
0246             end
0247         end
0248         if ~javaUnpackSuccessful
0249             %Copy content to temporary file
0250             tmpName = tempname();
0251             tmpFile = [tmpName '.gz'];
0252             hFileTmp = fopen(tmpFile, 'wb');
0253             
0254             if hFileTmp <= 0
0255                 error('Could not open temporary file for GZIP!');
0256             end
0257             
0258             fwrite(hFileTmp, compressedByteArray, 'uint8');
0259             fclose(hFileTmp);
0260             
0261             %Unzip with gzip
0262             gunzip(tmpFile);
0263             
0264             %Read the uncompressed file
0265             hFileTmp = fopen(tmpName, 'rb');
0266             if hFileTmp <= 0
0267                 error('Could not open unpacked file!');
0268             end
0269             cleanTmpFile = onCleanup(@() fclose(hFileTmp));
0270 
0271             cube = fread(hFileTmp, prod(metadata.cubeDim), metadata.datatype);
0272         end
0273     case {'bz2','bzip2'}
0274         error('Sorry: bzip compression not yet supported!');
0275     otherwise 
0276         error(['Undefined NRRD encoding scheme: ' nrrdMetaData.encoding]);
0277 end
0278 
0279 %maybe we need to correct the byte ordering (endian)
0280 if doSwapBytes
0281     swapbytes(cube);
0282 end
0283 
0284 if numel(metadata.cubeDim) > 1
0285     
0286     %first we shape the data into a cube
0287     cube = reshape(cube,metadata.cubeDim);
0288     
0289     %now we have to do the permutations, 2 1 3 ... is the MATLAB default
0290     %We create a transform matrix that transforms a permutation to MATLAB
0291     %default
0292     permutationTransformMatrix = diag(ones(metadata.dimension,1));
0293     permutationTransformMatrix(1:2,1:2) = flip(permutationTransformMatrix(1:2,1:2));
0294     
0295     applyPermutation = permutationTransformMatrix*metadata.axisPermutation';
0296     
0297     cube = permute(cube,applyPermutation);
0298 end
0299 
0300 
0301 
0302 end
0303 
0304 %This function wraps datatype definitions to the one used by MATLAB
0305 function datatype = getMATLABdataType(typestring)
0306 %The typedefinitions are taken directly from Section 5 of the NRRD
0307 %definition
0308 switch typestring
0309     case {'signed char','int8','int8_t'}
0310         datatype = 'int8';     
0311     case {'uchar','unsigned char','uint8','uint8_t'}
0312         datatype = 'uint8';        
0313     case {'short','short int','signed short','signed short int','int16','int16_t'}
0314         datatype = 'int16';        
0315     case {'ushort','unsigned short','unsigned short int','uint16','uint16_t'}
0316         datatype = 'uint16';        
0317     case {'int','signed int','int32','int32_t'}
0318         datatype = 'int32';        
0319     case {'uint','unsigned int','uint32','uint32_t'}
0320         datatype = 'uint32';        
0321     case {'longlong','long long','long long int','signed long long','signed long long int','int64','int64_t'}
0322         datatype = 'int64';        
0323     case {'ulonglong', 'unsigned long long', 'unsigned long long int','uint64','uint64_t'}
0324         datatype = 'uint64';        
0325     case 'float'
0326         datatype = 'single';        
0327     case 'double'
0328         datatype = 'double';        
0329     otherwise
0330         error('Could not identify datatype for NRRD data');
0331 end
0332 
0333 end
0334
matRad_readNRRD

Purpose

Synopsis

Description

Cross-reference information

Subfunctions

Source code
