createvmr.m

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

                spbs = sqrt(sum((slnc - sl1c) .^ 2)) / (DimZ - 1);
                slti = find(strcmp(p.MatrixHeaders, 'slice_thickness'));
                if isempty(slti)
                    error('UNSUPPORTED_PARREC');
                end
                pxz  = p.MatrixValues(1, slti(1));
                slgi = find(strcmp(p.MatrixHeaders, 'slice_gap'));
                if isempty(slgi)
                    error('UNSUPPORTED_PARREC');
                end
                gaps = p.MatrixValues(1, slgi(1));
                
                % possibly recalc gap?
                if abs(spbs - (pxz + gaps)) > 5e-4
                    gaps = spbs - pxz;
                end
                
                % initialize array size
                vdt(DimX, DimY, DimZ) = vdt(1);

                % get data access struct
                switch lower(p.RECDataType)
                    case {'dyn.slice'}
                        slio = p.RECData.Dyn(1).Slice;
                    case {'slice'}
                        slio = p.RECData.Slice;
                    otherwise
                        error('UNSUPPORTED_PARREC');
                end
                
                % iterate over slices
                for sc = 1:DimZ
                    vdt(:, :, sc) = slio(sc).IO(:, :);
                end
            catch
                error( ...
                    'BVQXtools:ErrorReadingFile', ...
                    'Error reading specified PAR/REC file.' ...
                );
            end
            
        % unsupported
        otherwise
            error( ...
                'BVQXtools:BadArgument', ...
                'Single file of this type not supported.' ...
            );
    end
else
    
    % also check extension
    switch lower(imafiles{1}(end-3:end))
        
        % DICOM / IMA
        case {'.dcm', '.ima'}
            try
                
                % assume number of files is number of slices
                DimZ = numel(imafiles);
                
                % read in first and last file
                im1 = BVQXfile(imafiles{1});
                cll{1} = im1;
                iml = BVQXfile(imafiles{end});
                cll{2} = iml;
                
                % check files
                if ~isBVQXfile(im1, 'dcm') || ...
                   ~isBVQXfile(iml, 'dcm') || ...
                    im1.DataLittleEndian ~= iml.DataLittleEndian
                    error('BAD_FILES');
                end
                
                % check some crucial settings
                dle = im1.DataLittleEndian;
                DimX = im1.Value('Columns');
                DimY = im1.Value('Rows');
                if ischar(DimX) || ...
                    ischar(DimY)
                    error('BAD_VALUES');
                end
                
                % check pixel data of first image
                DimS = DimX * DimY;
                IInt = im1.PixelData;
                if (~isa(IInt, 'uint16') && ...
                    ~isa(IInt, 'uint8')) || ...
                    numel(IInt) ~= DimS
                    error('WRONG_DIM_OR_TYPE');
                end
                
                % try to get some crucial parameters from the images
                try
                    iop = reshape(im1.Value('ImageOrientationPatient'), [3, 2])';
                    ipp = im1.Value('ImagePositionPatient');
                    ipl = iml.Value('ImagePositionPatient');
                    pxs = im1.Value('PixelSpacing');
                    pxz = im1.Value('SliceThickness');
                catch
                    iop = [0, 1, 0; 0, 0, -1];
                    ipp = [-(DimZ - 1) / 2, -DimX / 2, DimY / 2];
                end
                iop(3, :) = -cross(iop(1, :), iop(2, :));
                iop = iop ./ repmat(sqrt(sum(iop .* iop, 2)), [1, 3]);
                
                % compute BVQX settings
                sl1c = ipp + pxs(1) * (DimX / 2) * iop(1, :) + pxs(2) * (DimY / 2) * iop(2, :);
                slnc = ipl + pxs(1) * (DimX / 2) * iop(1, :) + pxs(2) * (DimY / 2) * iop(2, :);
                spbs = sqrt(sum((slnc - sl1c) .^ 2)) / (DimZ - 1);
                gaps = pxz - spbs;
                if abs(gaps) < 0.005
                    gaps = 0;
                end
                iop(4, :) = (slnc - sl1c) / (DimZ - 1);
                
                % initialize array size
                vdt(DimX, DimY, DimZ) = vdt(1);

                % iterate over slices
                for sc = 1:DimZ
                    fid = 0;
                    if dle
                        fid = fopen(imafiles{sc}, 'rb', 'ieee-le');
                    else
                        fid = fopen(imafiles{sc}, 'rb', 'ieee-be');
                    end
                    if fid < 3
                        error('FILE_OPEN_ERROR');
                    end
                    fseek(fid, -(12 + 2 * DimS), 1);
                    dtag = fread(fid, [1, 2], 'uint16=>double');
                    fread(fid, [1, 2], '*uint8');
                    dtls = fread(fid, [1, 1], 'uint16=>double');
                    dtll = fread(fid, [1, 1], 'uint32=>double');
                    if ~all(dtag == [32736, 16]) || ...
                        dtls ~= 0 || ...
                        dtll ~= (2 * DimS)
                        error('BAD_TAG_OR_LENGTH');
                    end
                    fcnt = fread(fid, [1, Inf], '*uint16');
                    if numel(fcnt) ~= DimS
                        error('FILE_TOO_SHORT');
                    end
                    vdt(:, :, sc) = reshape(fcnt, [DimX, DimY]);
                    fclose(fid);
                end

            catch
                clearbvqxobjects(cll);
                if fid > 0
                    fclose(fid);
                end
                error( ...
                    'BVQXtools:BVQXfileError', ...
                    'Error creating VMR from DICOM files.' ...
                );
            end
            
        % unsupported
        otherwise
            error( ...
                'BVQXtools:BadArgument', ...
                'Single file of this type not supported.' ...
            );
    end
end

% clear possible objects
clearbvqxobjects(cll);

% check orientation
iopr = corrcoef(iop(3, :), iop(4, :));
iopr = iopr(1, 2);
if abs(iopr) < 0.75
    error( ...
        'BVQXtools:InternalError', ...
        'Error getting image orientation.' ...
    );
end
if iopr < 0
    cnvt = 0;
else
    cnvt = 1;
end

% create object
vmr = bless(BVQXfile('new:vmr'), 1);

% set object properties
updstate = BVQXfile(0, 'updatestate');
BVQXfile(0, 'updatedisable', 'vmr');
vmr.DimX = DimX;
vmr.DimY = DimY;
vmr.DimZ = DimZ;
vmr.VMRData = [];
vmr.VMRData16 = [];
vmr.FramingCube = max([DimX, DimY, DimZ]);
fc = vmr.FramingCube / 2;
vmr.PosInfoVerified = 1;
vmr.OffsetX = max(0, ceil(fc - ceil(DimX / 2)));
vmr.OffsetY = max(0, ceil(fc - ceil(DimY / 2)));
vmr.OffsetZ = max(0, ceil(fc - ceil(DimZ / 2)));
vmr.Slice1CenterX = sl1c(1);
vmr.Slice1CenterY = sl1c(2);
vmr.Slice1CenterZ = sl1c(3);
vmr.SliceNCenterX = slnc(1);
vmr.SliceNCenterY = slnc(2);
vmr.SliceNCenterZ = slnc(3);
vmr.RowDirX = iop(1, 1);
vmr.RowDirY = iop(1, 2);
vmr.RowDirZ = iop(1, 3);
vmr.ColDirX = iop(2, 1);
vmr.ColDirY = iop(2, 2);
vmr.ColDirZ = iop(2, 3);
vmr.NRows = DimX;
vmr.NCols = DimY;
vmr.FoVRows = DimX * pxs(1);
vmr.FoVCols = DimY * pxs(2);
vmr.SliceThickness = pxz;
vmr.GapThickness = gaps;
vmr.Convention = cnvt;
vmr.VoxResX = pxs(1);
vmr.VoxResY = pxs(2);
vmr.VoxResZ = pxz + gaps;
vmr.VoxResInTalairach = intal;
vmr.VoxResVerified = true;

% set data
switch (lower(class(vdt)))
    case {'int8', 'int16'}
        vdt = uint16(single(vdt) - single(min(vdt(:))));
    case {'uint8'}
        vdt = uint16(vdt);
    case {'uint32'}
        if max(vdt(:)) > 32767
            vdt = uint16(32767 * (single(vdt) / max(vdt(:))));
        else
            vdt = uint16(vdt);
        end
    case {'int32', 'single', 'double'}
        vdt = single(vdt) - single(min(vdt(:)));
        if max(vdt(:)) > 32767
            vdt = uint16(32767 * (vdt / max(vdt(:))));
        else
            vdt = uint16(vdt);
        end
    case {'uint16'}
        % do nothing
    otherwise
        error( ...
            'BVQXtools:BadFileContent', ...
            'Invalid VMR creation datatype.' ...
        );
end
vmr.VMRData16 = vdt;
vmr.MinOriginalValue = double(min(vdt(:)));
vmr.MeanOriginalValue = round(double(mean(vdt(:))));
vmr.MaxOriginalValue = double(max(vdt(:)));

% try to perform good limitting
[hn, hx] = hist(single(vdt(:)), single(vmr.MinOriginalValue:vmr.MaxOriginalValue));
hx = uint16(hx);
cn = cumsum(hn);
lt = find(cn > 0.01 * cn(end));
ut = find(cn > 0.99 * cn(end));
if isempty(lt)
    lt = 1;
else
    lt = lt(1);
end
if isempty(ut)
    ut = numel(cn);
else
    ut = ut(1);
end
vdt(vdt < hx(lt)) = hx(lt);
vdt(vdt > hx(ut)) = hx(ut);
vdt = single(vdt) - single(hx(lt));
vmr.VMRData = uint8(225 * (vdt ./ max(vdt(:))));
BVQXfile(0, 'updatestate', updstate);