利用靶区勾画的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
上一篇: Mathematica斐波那契线搜索代码