DICOM dimensions in Matlab matrix (all frames end in the last dimension of the array)

, . , , , .

:

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);

:

%% 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]);

Ive , mathworks.