DICOM dimensions in matlab array (all frames end up in last dimension of array)

Question

In one of my GUIs I load DICOM images. Sometimes they are just a volume and another dimension and when I load them in Matlab everything ends up where I want it.

handles.inf = dicominfo([filepath filename]);
handles.dat = dicomread(handles.inf);
size(handles.dat)

ans = 128 128 128 512

For an 128 by 128 by 128 volume at 512 timepoints for example (actually the third dimension would not even be 128, the third dimension is stacks, of which I don't know what it is). However sometimes There are more dimensions in the dicom, but the reader just puts all of them in the fourth dimension.

handles.inf = dicominfo([filepath filename]);
handles.dat = dicomread(handles.inf);
size(handles.dat)

ans = 128 128 1  4082

For a single 128 by 128 slice with 512 timepoints, two echoes and magnitude, phase, real and imaginary data for example.

It is then very hard to unscramble them. Manually I can do this for every DICOM I load but when in a GUI I would want to have a general approach that just creates a dimension in the array for each dimension in the dicom.

This is especially important not just for data-analysis, but also to transform the coordinates from image space to patient space. My own approach was to look at the header, but there's no guarantee that certain entries will work, and the order in which they are applied I can't find. The header entries I found so far:

inf.Rows;%inf.width;%inf.MRAcquisitionFrequencyEncodingSteps;%inf.MRAcquisitionPhaseEncodingStepsInPlane
inf.Columns;% inf.height; % inf.NumberOfKSpaceTrajectories;
inf.MRSeriesNrOfSlices
inf.MRSeriesNrOfEchoes
inf.MRSeriesNrOfDynamicScans
inf.MRSeriesNrOfPhases
inf.MRSeriesReconstructionNumber % not sure about this one
inf.MRSeriesNrOfDiffBValues
inf.MRSeriesNrOfDiffGradOrients
inf.MRSeriesNrOfLabelTypes

reshapeddat = reshape(dat, [all dimension sizes from header here]);

I'm not sure how to check if I've got all variables and what the right order for the reshape. Anybody knows of a sure-fire way to get all dimension sizes from the DICOM header and the order in which they are stacked?

Answer 1

Ok so I now manually go by all possible dimensions. When a stack also contains reconstructed data which has less dimensions than the rest, remove those first.

This is how I check the dimensions:

info = dicominfo(filename);
datorig = dicomread(filename);

%dimension sizes per frame
nrX = double(info.Rows);              %similar nX;% info.width;% info.MRAcquisitionFrequencyEncodingSteps;% info.MRAcquisitionPhaseEncodingStepsInPlane
nrY = double(info.Columns);           %similar nY;% info.height;% info.NumberOfKSpaceTrajectories;

%dimensions between frames
nrEcho = double(info.MRSeriesNrOfEchoes);
nrDyn = double(info.MRSeriesNrOfDynamicScans);
nrPhase = double(info.MRSeriesNrOfPhases);
nrSlice = double(info.MRSeriesNrOfSlices);           %no per frame struct entry, calculate from offset.

%nr of frames
nrFrame = double(info.NumberOfFrames);

nrSeq = 1;                                   % nSeq not sure how to interpret this, wheres the per frame struct entry?
nrBval = double(info.MRSeriesNrOfDiffBValues);        % nB
nrGrad = double(info.MRSeriesNrOfDiffGradOrients);    % info.MRSeriesNrOfPhaseContrastDirctns;%similar nGrad?
nrASL = 1;                                   % info.MRSeriesNrOfLabelTypes;%per frame struct entry?

imtype = cell(1, nrFrame);
for ii = 1:nrFrame
    %imtype(ii) = eval(sprintf('info.PerFrameFunctionalGroupsSequence.Item_%d.PrivatePerFrameSq.Item_1.MRImageTypeMR', ii));
    imtype{ii} = num2str(eval(sprintf('info.PerFrameFunctionalGroupsSequence.Item_%d.PrivatePerFrameSq.Item_1.MRImageTypeMR', ii)));
end

imType = unique(imtype, 'stable');
nrType = length(imType);

This is how I reformat the dimensions:

%% count length of same dimension positions from start
if nrEcho > 1
    for ii = 1:nrFrame
        imecno(ii) = eval(sprintf('inf.PerFrameFunctionalGroupsSequence.Item_%d.PrivatePerFrameSq.Item_1.EchoNumber', ii));
    end
    lenEcho = find(imecno ~= imecno(1), 1, 'first') - 1;
else
    lenEcho = nrFrame;
end

if nrDyn > 1
    for ii = 1:nrFrame
        imdynno(ii) = eval(sprintf('inf.PerFrameFunctionalGroupsSequence.Item_%d.PrivatePerFrameSq.Item_1.TemporalPositionIdentifier', ii));
    end
    lenDyn = find(imdynno ~= imdynno(1), 1, 'first') - 1;
else
    lenDyn = nrFrame;
end

if nrPhase > 1
    for ii = 1:nrFrame
        imphno(ii) = eval(sprintf('inf.PerFrameFunctionalGroupsSequence.Item_%d.PrivatePerFrameSq.Item_1.MRImagePhaseNumber', ii));
    end
    lenPhase = find(imphno~=imphno(1), 1, 'first') - 1;
else
    lenPhase = nrFrame;
end


if nrType > 1
    q = 1;
    imtyno(1, 1) = q;
    for ii = 2:nrFrame        
        if imtype{:, ii-1} ~= imtype{:, (ii)}
            q = q+1;
        end
        imtyno(1, ii) = q;
        %for jj = 1:nrType
            %if imtype{:,ii} == imType{:,jj} 
            %    imtyno(1, ii) = jj;
            %end
        %end
    end
    if q ~= nrType
        nrType = q;
    end
    lenType = find(imtyno ~= imtyno(1), 1, 'first') - 1;
else
    lenType = nrFrame;
end

% slices are not numbered per frame, so get this indirectly from location
% currently not sure if working for all angulations
for ii = 1:nrFrame
    imslice(:,ii) =  -eval(['inf.PerFrameFunctionalGroupsSequence.Item_',sprintf('%d', ii),'.PlanePositionSequence.Item_1.ImagePositionPatient']);
end

% stdsl = std(imslice,[],2); --> Assumption
% dirsl = max(find(stdsl == max(stdsl)));
imslices = unique(imslice', 'rows')';
if nrSlice > 1
    for ii = 1:nrFrame
        for jj = 1:nrSlice
            if imslice(:,ii) == imslices(:,nrSlice - (jj-1)); %dirsl or :?
                imslno(1, ii) = jj;
            end
        end
    end
    lenSlice = find(imslno~=imslno(1), 1, 'first')-1;
else
    lenSlice = nrFrame;
end

if nrBval > 1
    for ii = 1:nrFrame
        imbno(ii) = eval(sprintf('inf.PerFrameFunctionalGroupsSequence.Item_%d.PrivatePerFrameSq.Item_1.MRImageDiffBValueNumber', ii));
    end
    lenBval = find(imbno~=imbno(1), 1, 'first') - 1;
else
    lenBval = nrFrame;
end

if nrGrad > 1
    for ii = 1:nrFrame
        imgradno(ii) = eval(sprintf('inf.PerFrameFunctionalGroupsSequence.Item_%d.PrivatePerFrameSq.Item_1.MRImageGradientOrientationNumber', ii));
    end
    lenGrad = find(imgradno~=imgradno(1), 1, 'first')-1;
else
    lenGrad = inf.NumberOfFrames;
end

lenSeq = nrFrame; % dont know how to get this information per frame, in my case always one
lenASL = nrFrame; % dont know how to get this information per frame, in my case always one

%this would have been the goal format
goaldim = [nrSlice nrEcho nrDyn nrPhase nrType nrSeq nrBval nrGrad nrASL];             % what we want
goallen = [lenSlice lenEcho lenDyn lenPhase lenType lenSeq lenBval lenGrad lenASL];    % what they are

[~, permIX] = sort(goallen);

dicomdim = zeros(1, 9);
for ii = 1:9
    dicomdim(1, ii) = goaldim(permIX(ii));
end

dicomdim = [nrX nrY dicomdim];
%for any possible zero dimensions from header use a 1 instead
dicomdim(find(dicomdim == 0)) = 1;

newdat = reshape(dat, dicomdim);

newdim = size(newdat);
newnonzero = length(newdim(3:end));
goalnonzero = permIX(1:newnonzero);
[dummyy, goalIX] = sort(goalnonzero);
goalIX = [1 2 goalIX+2]; 
newdat = permute(newdat, goalIX);
newdat = reshape(newdat, [nrX nrY goaldim]);

When Ive used this as a function for a longer period and debugged it a bit I might post in on the file exchange of mathworks.