210812-Matlab读取StrctrSets.dcm文件(轮廓dcm，210918更新）

引文：记录Matlab读取StrctrSets.dcm文件（轮廓dcm)，但这里读取的是坐标信息，还需要与原始CT dicom结合，将坐标转化为矩阵，这里没有尝试了。（210918更新）

参考：How to display DICOM RT structure

clc;clear; close all;
info=dicominfo('Y201040_StrctrSets.dcm');
% contour = dicomContours(info);  % R2020b新版中有，2016b无
% contour.ROIs
structures=fieldnames(info.ROIContourSequence); %Lists all the structures
num_cont=size(structures,1);                        %Says how many it is
item_pos_all=[];
for n=1:num_cont  %Loops through all structures    
    curr_struct=structures{n};    
    %Lists all the ROI-layers within a structure
    ROIlayers=fieldnames(info.ROIContourSequence.(curr_struct).ContourSequence);
    num_layers=size(ROIlayers,1);        
    for i=1:num_layers     
	    currROIlayer=ROIlayers{i};
	    item_pos=info.ROIContourSequence.(curr_struct).ContourSequence.(currROIlayer).ContourData;   % 读取每层轮廓所有连接点的坐标
	    item_pos=reshape(item_pos, [3 length(item_pos)/3]);
	    item_pos_all=[item_pos_all, item_pos];    
    end       
    figure(n)
    hold on
    plot3(item_pos_all(1,:), item_pos_all(2,:), item_pos_all(3,:)) 
end

后续210812

后来在上面基础上，自己将坐标转化为矩阵，代码如下：

function [masks,masks_names]=getStructContour(patient_path)
%getStructContour Get structure contour and names from patient_path
%   patient_path: need only contains sorted CT dicom and structure
%   dicom(structure dicom is listed at last)
    dicom_names=dir(patient_path);  % List all dicom names
    struct_name=dicom_names(end).name;  % Get structure dicom name
    struct_path=fullfile(patient_path,struct_name);
    benchmark_dicom_name=dicom_names(end-1).name;  % Get benchmark dicom name
    benchmark_dicom_path=fullfile(patient_path,benchmark_dicom_name);
    
    info_struct=dicominfo(struct_path);  % Get structure dicom info
    structures=fieldnames(info_struct.ROIContourSequence); %Lists all the structures
    num_struct=size(structures,1);                        %Says how many it is
    info_dcm=dicominfo(benchmark_dicom_path);  % Get benchmark dicom info
    masks=zeros(info_dcm.Rows,info_dcm.Columns,length(dir(patient_path))-3,num_struct);  % 4D mask
    masks_names=cell(num_struct,1);  % store contour(here alse name mask) name by cell
    
    for ii=1:num_struct  % Loops through all structures    
        curr_struct=structures{ii}; 
        disp([patient_path,'---',info_struct.StructureSetROISequence.(curr_struct).ROIName]);
        masks_names{ii}=[info_struct.StructureSetROISequence.(curr_struct).ROIName];  % store masks_name
        % Lists all the ROI-layers within a structure
        ROIlayers=fieldnames(info_struct.ROIContourSequence.(curr_struct).ContourSequence);  
        num_slice_bystruct=size(ROIlayers,1);
        if  num_slice_bystruct~= size(masks,3)  % check slice num by struct == num by dicom num
            fprintf('%s    %s   number ~= dicom number!\n',patient_path,masks_names{ii})
%             return
        end
        for jj=1:num_slice_bystruct %size(masks,3)  % Loop all slices in one structure 
            currROIlayer=ROIlayers{jj}; 
            % Get x,y,z points coordinates in one slice
            index_pos=info_struct.ROIContourSequence.(curr_struct).ContourSequence.(currROIlayer).ContourData;   
            index_pos=reshape(index_pos, [3 length(index_pos)/3]);
            index_pos=bsxfun(@minus, index_pos, info_dcm.ImagePositionPatient);
            index_pos=bsxfun(@rdivide, index_pos, [info_dcm.PixelSpacing;info_dcm.SliceThickness]); % from points to matrix
            index_pos=index_pos+1;  % chang index start from 1
            index_pos(index_pos<1)=1;  % prevent drawing mask out of boundaries
            index_pos=single(round(index_pos));  % round off
            
            mask=zeros(info_dcm.Rows,info_dcm.Columns);
            num_points=size(index_pos,2);
            
            for kk=1:num_points
                p1 = index_pos(1:2,kk);                      
                if kk~=num_points
                    p2 = index_pos(1:2,kk+1);
                else
                    p2 = index_pos(1:2,1);  % connect the head and tail
                end
                num=max(abs(p1(1)-p2(1)), abs(p1(2)-p2(2)))*2;  % use the big number of x and y(understand by slope of 1)
                x=int16(linspace(p1(1),p2(1),num));  % use linspace instead of line equation, Grea
                y=int16(linspace(p1(2),p2(2),num));
                mask(sub2ind(size(mask), y, x)) = 1;  % mask(y,x) is wrong!          
            end
%             imshow(mask,[])
            masks(:,:,index_pos(3,1),ii)=masks(:,:,index_pos(3,1),ii) | mask;   % this is the simplest method                
            
        end
    end
    ind=size(masks,3):-1:1;  % resort masks to fit with dicom images
    masks=masks(:,:,ind,:);
    % for ii=1:num_cont
    %     for jj=1:size(masks,3)
    %         imshow(masks(:,:,jj,ii),[]), title(masks_names(ii),'Interpreter','none')
    %         pause(0.1)
    %     end
    % end
end

210918更新

这个网站https://github.com/LiJiongliang/dicom-matlab-tools提供了图像、剂量、轮廓的读取，自己在此基础上自己研读，还是收获颇多。这里可以参照原文，注意坐标的读取与选择。