利用靶区勾画的RT struct文件分析相关联的Dicom图像中肿瘤部分的CT(HU)值分布

利用已经进行靶区勾画的RT struct文件分析相关联的Dicom图像中肿瘤部分的CT(HU)值分布。

% =================================================================
% 根据RT struct的靶区勾画区域分析相应dicom图像中肿瘤位置的CT HU值的分布
% =================================================================

%% 读取RT struct 靶区头文件
% rt struct文件路径
RSInfo = dicominfo('E:\RTFolder\rt.dcm');
% dicom系列文件所在文件夹
dcmPath = 'E:\ DicomFolder\';

numberOfContours = size(fieldnames(RSInfo.ROIContourSequence.Item_1.ContourSequence),1); % 该靶区分布在多少层 %已确定标注中Item_1为肿瘤区域勾画
maskData = zeros(512,512,numberOfContours);
tumorFile = maskData;

%% 遍历每一个切片上的勾画区域
for k=1:numberOfContours
    rfContent = RSInfo.ROIContourSequence.Item_1.ContourSequence.(['Item_' num2str(k)]);
    
    %% 读取该靶区所在的CT切片的信息
    dcmName = rfContent.ContourImageSequence.Item_1.ReferencedSOPInstanceUID;
    DCMInfo = dicominfo(strcat(dcmPath, dcmName, '.dcm')); 
    dcmFile = dicomread(strcat(dcmPath, dcmName, '.dcm'));
    dcmHU = dcmFile.* DCMInfo.RescaleSlope + DCMInfo.RescaleIntercept; % 将C灰度值转换为HU值
    dcmOrigin = DCMInfo.ImagePositionPatient; % 网格原点在世界坐标系的位置
    dcmSpacing = DCMInfo.PixelSpacing; %采样间隔
 
    
    numberOfPoints = rfContent.NumberOfContourPoints; % 该层靶区曲线点数
    conData = zeros(numberOfPoints,3); % 存储靶区曲线各点的世界坐标
    pointData = zeros(numberOfPoints,2); % 存储靶区曲线各点的网格体素坐标
    
    %% 将靶区勾画的曲线坐标由世界坐标系转换为网格体素坐标
    for jj = 1:numberOfPoints
        ii = (jj-1)*3 ;
        conData(jj,1) = rfContent.ContourData(ii+1,1); %轮廓世界坐标系
        conData(jj,2) = rfContent.ContourData(ii+2,1);
        conData(jj,3) = rfContent.ContourData(ii+3,1);       
        pointData(jj,1) = round( (conData(jj,1) - dcmOrigin(1,1))/dcmSpacing(1,1) ); %轮廓X坐标
        pointData(jj,2) = round( (conData(jj,2) - dcmOrigin(2,1))/dcmSpacing(2,1) ); %轮廓Y坐标   
    end
     
   %% 判断每一个切面上所有的点是否在曲线内部，在maskData相应位置=1,不在=0； 
    x = zeros(1,512);
    y = x;
    for i = 1:512
        x(i,:) = i;
        y(i,:) = 1:512;   
    end
    in = inpolygon(x,y,pointData(:,1)', pointData(:,2)');
    maskData(:,:,k) = in;
    
    tumorFile(:,:,k) = dcmHU;
    tumorFile(~maskData) = 0;  % 只有肿瘤部分的原始图像
    
    %% 显示原图像和肿瘤部分
    figure;
    imshow(int8(dcmHU)); % 显示整体图像 %将int16转化为int8，显示出来易于观察
    hold on;
    plot(pointData(:,1),pointData(:,2),x(in),y(in),'.r'); % 显示肿瘤
        
end

%% 分析肿瘤部分的HU值的频数和频率
HUFqc = tabulate(tumorFile(:));
% 去掉频数低于某一阈值的统计数据
HUFqcNew = zeros(100,size(HUFqc,2),size(HUFqc,3));
rowNum = 0;
for j = 1:size(HUFqc,1)
    if HUFqc(j,2) > 100
        rowNum = rowNum + 1;
        HUFqcNew(rowNum,:) = HUFqc(j,:);
    end
    
end
% 绘制直方图
% tumorHist = tumorFile(:);
% figure;
% hist(tumorHist,10); %绘制分布直方图

参考： 

http://blog.csdn.net/sunyao_123/article/details/73655957  sunyao_123的《DICOM靶区头文件解析》

http://blog.csdn.net/u013635029/article/details/72957944?locationNum=2&fps=1  Hade_H的《常见医疗扫描图像处理步骤》

如有任何不正确的地方，欢迎指正。