欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

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

程序员文章站 2022-03-02 10:41:54
...

利用靶区勾画的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







相关标签: matlab CT影像