dcmio.m

toolbox of BVQX, This is the access between BV and matlab. It will help you to analysis data from BV

function [varargout] = dcmio(varargin)
% dcmio  - read DICOM files
%
% FORMAT:       dcmobject = dcmio(filename)
%
% Input fields:
%
%       filename    filename of DICOM file to read
%
% Output fields:
%
%       dcmobject   BVQXfile object
%
% See also BVQXfile

% Version:  v0.6d
% Build:    7020711
% Date:     Feb-07 2007, 11:45 AM CET
% Author:   Jochen Weber, Brain Innovation, B.V., Maastricht, NL
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% persistent VR dict and empty file
persistent my_dicom_vr my_dicom_dic my_newdcm;
if isempty(my_dicom_vr)
    my_dicom_vr = dicom_vr;
    my_dicom_dic = dicom_dic(4);
    my_newdcm = BVQXfile(0, 'newcont', 'dcm');
end
vrSQ = my_dicom_vr.SQ;
vrUN = my_dicom_vr.UN;

% argument check
if nargin < 1 || ...
    ~ischar(varargin{1}) || ...
    isempty(varargin{1})
    error( ...
        'BVQXtools:BadArgument', ...
        'Bad or missing argument for dcmio.' ...
    );
end
filename = varargin{1}(:)';

% default options
writemode = false;

% dcm content structure given
if nargin > 1 && ...
    isBVQXfile(varargin{2}, 'dcm') && ...
    numel(varargin{2}) == 1
    try
        dcmcont = getcont(varargin{2});
    catch
        error( ...
            'BVQXtools:BVQXfileError', ...
            'DICOM object no longer available.' ...
        );
    end
    writemode = true;
end

% reading DCM content
if ~writemode

    % create new object
    dcmcont = my_newdcm;

    % check file (opening in little endian syntax by default)
    fid = fopen(filename, 'rb', 'ieee-le');
    if fid < 1
        error( ...
            'BVQXtools:FileNotReadable', ...
            'File not readable: ''%s''.', ...
            filename ...
        );
    end

    % get file size
    fseek(fid, 0, 1);
    flen = ftell(fid);
    flenmin = flen - 8;
    fseek(fid, 0, -1);
    if flen < 132
        fclose(fid);
        error( ...
            'BVQXtools:InvalidFilesize', ...
            'File too short: ''%s''.', ...
            filename ...
        );
    end

    % get file content (little endian)
    fc_uint8 = fread(fid, [1, Inf], '*uint8');
    fseek(fid, 0, -1);
    fc_uint16_le = fread(fid, [1, Inf], '*uint16');
    fseek(fid, 0, -1);

    % get file content (big endian)
    fidbe = fopen(filename, 'rb', 'ieee-be');
    fc_uint16_be = fread(fidbe, [1, Inf], '*uint16');
    fseek(fidbe, 0, -1);

% writing BFF content
else

    % check open writable file
    if dcmcont.MetaLittleEndian
        fid = fopen(filename, 'wb', 'ieee-le');
    else
        fid = fopen(filename, 'wb', 'ieee-be');
    end
    if fid < 1
        error( ...
            'BVQXtools:FileNotWritable', ...
            'File not writable: ''%s''.', ...
            filename ...
        );
    end
    fclose(fid);
    
    % however, still not yet implemented
    error( ...
        'BVQXtools:NotYetImplemented', ...
        'Writing of DICOM files not yet implemented.' ...
    );

end

% counters
MetaKeyCount = 1;
DataKeyCount = 0;

% try used for forcemode (on indentation !!!)
try % forcemode

    % reading
    if ~writemode

        % detect magic token
        dcmcont.Magic = 'DICM';
        dcmcont.Preamble = true;

        % if magic token present
        if all(fc_uint8(129:132) == [68, 73, 67, 77])

            % start reading after token
            fpos = 133;

        % otherwise
        else

            % reset defaults and start at char 1!
            dcmcont.Magic = '';
            dcmcont.Preamble = false;
            fpos = 1;
        end

        % defaults, detection
        dcmcont.MetaTSExplicit = false;
        dcmcont.MetaLittleEndian = true;

        % if Magic token found
        if dcmcont.Preamble

            % initialize meta struct
            MetaStr = dcmcont.Meta;

            % check implicit explicit
            vrcheck = fc_uint8(fpos+4:fpos+5);
            dcmcont.MetaTSExplicit = all(vrcheck > 64 & vrcheck < 91);

            % check endian only on explicit VR syntax
            MetaExplicit = dcmcont.MetaTSExplicit;
            if MetaExplicit
                lecheck = fc_uint8(fpos:fpos+4);
                if lecheck(1) < lecheck(2)
                    dcmcont.MetaLittleEndian = false;
                elseif lecheck(1) == lecheck(2)
                    if lecheck(3) < lecheck(4)
                        dcmcont.MetaLittleEndian = false;
                    elseif lecheck(3) == lecheck(4)
                        error( ...
                            'BVQXtools:FileDetectionFailed', ...
                            'Error detecting Endian type in Meta header.' ...
                        );
                    end
                end
            end
            MetaLEndian = dcmcont.MetaLittleEndian;
            if MetaLEndian
                rfid = fid;
            else
                rfid = fidbe;
            end

            % read until key is over group 7
            while true

                % get correct key format
                if MetaLEndian
                    DicomKey = double(fc_uint16_le(round(fpos/2):round(fpos/2)+1));
                else
                    DicomKey = double(fc_uint16_be(round(fpos/2):round(fpos/2)+1));
                end

                % break on condition
                if DicomKey(1) > 7
                    break;
                end
                fpos = fpos + 4;
                MetaStr(MetaKeyCount).Key = DicomKey;
                if DicomKey(1) > 65533
                    MetaSequence = true;
                    DataVR = vrSQ;
                else
                    MetaSequence = false;
                end

                % VR given ?
                if ~MetaSequence && ...
                    MetaExplicit
                    DataVR = upper(char(fc_uint8(fpos:fpos+1)));
                    fpos = fpos + 2;
                    if ~isfield(my_dicom_vr, DataVR)
                        error( ...
                            'BVQXtools:InvalidToken', ...
                            'Invalid VR found: %s.', ...
                            DataVR ...
                        );
                    end
                    DataVR = my_dicom_vr.(DataVR);
                elseif ~MetaSequence
                    DataVR = vrUN;
                end
                MetaStr(MetaKeyCount).VR = DataVR.tag;
                if DataVR.length(2) > 32767
                    DataShortVLength = false;
                else
                    DataShortVLength = true;
                end

                % length
                if MetaLEndian
                    DataVLengthShort = double(fc_uint16_le(round(fpos/2)));
                else
                    DataVLengthShort = double(fc_uint16_be(round(fpos/2)));
                end
                if ~DataShortVLength && ...
                   ~MetaSequence && ...
                    MetaExplicit && ...
                    DataVLengthShort > 0
                    error( ...
                        'BVQXtools:InvalidFileContent', ...
                        'Invalid 16-bit VL for given VR.' ...
                    );
                elseif MetaSequence
                    fpos = fpos - 2;
                end
                if ~DataShortVLength
                    if MetaLEndian
                        DataVLengthLong = 65536 * ...
                            double(fc_uint16_le(round(fpos/2)+2)) + ...
                            double(fc_uint16_le(round(fpos/2)+1));
                    else
                        DataVLengthLong = 65536 * ...
                            double(fc_uint16_le(round(fpos/2)+1)) + ...
                            double(fc_uint16_le(round(fpos/2)+2));
                    end
                    fpos = fpos + 6;
                    dvl = DataVLengthLong;
                else
                    fpos = fpos + 2;
                    DataVLengthLong = NaN;
                    dvl = DataVLengthShort;
                end
                MetaStr(MetaKeyCount).VLShort = DataVLengthShort;
                MetaStr(MetaKeyCount).VLLong  = DataVLengthLong;
                dvl = 2 * round(dvl/2);

                % get value according to type
                ddt = lower(DataVR.datatype);
                if ~strcmp(ddt, 'sequence')

                    % go to and update position
                    fseek(rfid, fpos - 1, -1);
                    fpos = fpos + dvl;