matlab数组中的DICOM维度（所有帧都以数组的最后一维结束）

Question

在我的一个 GUI 中，我加载了 DICOM 图像。有时它们只是一个体积和另一个维度，当我将它们加载到 Matlab 中时，一切都会在我想要的地方结束。

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

ans = 128 128 128 512

例如，对于 512 个时间点的 128 x 128 x 128 卷（实际上第三维甚至不会是 128，第三维是堆栈，我不知道它是什么）。然而有时dicom中的维度更多，但读者只是将它们全部放在第四维度中。

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

ans = 128 128 1  4082

对于具有 512 个时间点、两个回波和幅度、相位、实部和虚部数据的单个 128 x 128 切片。

然后很难解读它们。手动我可以为我加载的每个 DICOM 执行此操作，但是在 GUI 中我希望有一种通用方法，只需在数组中为 dicom 中的每个维度创建一个维度。

这不仅对数据分析特别重要，而且对transform 从图像空间到患者空间的坐标也特别重要。我自己的方法是查看标题，但不能保证某些条目会起作用，而且我找不到它们的应用顺序。到目前为止我找到的标题条目：

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

我不确定如何检查我是否拥有所有变量以及重塑的正确顺序。有人知道从 DICOM 标头中获取所有尺寸大小及其堆叠顺序的可靠方法吗？

Answer 1

好的，我现在手动遍历所有可能的维度。当堆栈还包含维数小于其余数据的重构数据时，请先删除这些数据。

这是我检查尺寸的方式：

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

当我将它作为一个函数使用了很长时间并对其进行了一些调试时，我可能会在 mathworks 的文件交换中发帖。