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

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

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

%% 读取RT struct 靶区头文件

% 设置路径
dstPath = 'E:\dstFolder\';
srcPath = 'E:\srcFolder\';
patient_dir = dir([srcPath,'*']);
for patNum =3:length(patient_dir)  %遍历每个病人
% for patNum = 3:3
    image_dir = dir([srcPath, patient_dir(patNum).name, '\*']);
    for imageNum = 3:length(image_dir)  %遍历该病人的不同的扫描CT
        ct_dir = dir([image_dir(imageNum).folder,'\', image_dir(imageNum).name, '\*', '\*']);
        % dicom系列文件所在文件夹
        dcmPath = [ct_dir(3).folder, '\', ct_dir(3).name, '\'];
        % rt struct文件路径
        RSInfoPath = dir([ct_dir(4).folder, '\', ct_dir(4).name, '\*.dcm']);
        RSInfo = dicominfo(strcat(RSInfoPath.folder, '\', RSInfoPath.name));
        
        % 开始处理
        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) > 80
                rowNum = rowNum + 1;
                HUFqcNew(rowNum,:) = HUFqc(j,:);
            end

        end
        dst_dir = strcat(dstPath,'\', patient_dir(patNum).name);
        if ~exist(dst_dir)
            mkdir(dst_dir);
        end
        dst_name1 = strcat(dst_dir, '\',image_dir(imageNum).name, '_HUFqc.xls');
        dst_name2 = strcat(dst_dir, '\',image_dir(imageNum).name, '_HUFqcNew.xls');
%         if exist(dst_name1)
%             delete dst_name1;
%         end
%         if exist(dst_name2)
%             delete dst_name2;
%         end
        
        xlswrite(dst_name1, HUFqc);
        xlswrite(dst_name2,HUFqcNew);
        % 绘制直方图
        % tumorHist = tumorFile(:);
        % figure;
        % hist(tumorHist,10); %绘制分布直方图
        
        
        
        
    end
    
end