CT dicom 5mm层厚插值成1mm-20201110

放射科的CT大多数是5 mm,有时需要插值成1 mm,这里是医工所的曹老师写的代码,在此做个记录,学习一下。

clear;clc;
%程序功能:读取dicom文件,并将其插值为1mm
folderPath = 'Y180521\';
outFolderPath = 'out\';

folderInfo = dir(folderPath);
fileNum = size(folderInfo,1)-2;
outFileNum = 0;
fileIndex = [];
instanceNumber = [];
for i = 1:fileNum
    %读取相邻两张图片
    fileName = folderInfo(i+2).name;
    filePath = strcat(folderPath, fileName);
%     ImgB = dicomread(filePath);  % 这里似乎不需要,只需要头文件信息
    ImgInfoB = dicominfo(filePath);%存储信息
%     disp(ImgInfoB.InstanceNumber);
    instanceNumber(end+1,1) = ImgInfoB.InstanceNumber;  % end+1后的,1可省略
end
startInstanceNumber = min(instanceNumber);  % 细心,注意设置起始number
for i = 0:(fileNum-2)
    %读取相邻两张图片
    fileIndexB = find(instanceNumber == (startInstanceNumber + i));    
    fileNameB = folderInfo(fileIndexB+2).name;
    filePathB = strcat(folderPath, fileNameB);
    ImgB = dicomread(filePathB);
    ImgInfoB = dicominfo(filePathB);%存储信息
    
    fileIndexT = find(instanceNumber == (startInstanceNumber + i + 1));
    fileNameT = folderInfo(fileIndexT+2).name;
    filePathT = strcat(folderPath, fileNameT);
    ImgT = dicomread(filePathT);
    ImgInfoT = dicominfo(filePathT);%存储信息

    %线性插值
    ImgD = double(ImgT - ImgB) * 0.2;
    Img1 = int16(round(ImgD * 1 + double(ImgB)));
    Img2 = int16(round(ImgD * 2 + double(ImgB)));
    Img3 = int16(round(ImgD * 3 + double(ImgB)));
    Img4 = int16(round(ImgD * 4 + double(ImgB)));
    
    %保存结果
    filePath = strcat( outFolderPath,num2str(i*5),'.DCM');
    ImgInfoB.InstanceNumber = i * 5 + 1;
    ImgInfoB.SliceThickness = 1;
    ImgInfoB.Filename = filePath;
    dicomwrite(ImgB,filePath,ImgInfoB);%写入Dicom图像

    filePath = strcat( outFolderPath,num2str(i*5+1),'.DCM');
    ImgInfoB.InstanceNumber = i * 5 + 2;
    ImgInfoB.SliceThickness = 1;
    ImgInfoB.Filename = filePath;
    ImgInfoB.SliceLocation = ImgInfoB.SliceLocation - 1;
    ImgInfoB.ImagePositionPatient(3,1) = ImgInfoB.ImagePositionPatient(3,1)-1;
    dicomwrite(Img1,filePath,ImgInfoB);%写入Dicom图像
    
    filePath = strcat( outFolderPath,num2str(i*5+2),'.DCM');
    ImgInfoB.InstanceNumber = i * 5 + 3;
    ImgInfoB.SliceThickness = 1;
    ImgInfoB.Filename = filePath;
    ImgInfoB.SliceLocation = ImgInfoB.SliceLocation - 1;
    ImgInfoB.ImagePositionPatient(3,1) = ImgInfoB.ImagePositionPatient(3,1)-1;
    dicomwrite(Img2,filePath,ImgInfoB);%写入Dicom图像
    
    filePath = strcat( outFolderPath,num2str(i*5+3),'.DCM');
    ImgInfoB.InstanceNumber = i * 5 + 4;
    ImgInfoB.SliceThickness = 1;
    ImgInfoB.Filename = filePath;
    ImgInfoB.SliceLocation = ImgInfoB.SliceLocation - 1;
    ImgInfoB.ImagePositionPatient(3,1) = ImgInfoB.ImagePositionPatient(3,1)-1;
    dicomwrite(Img3,filePath,ImgInfoB);%写入Dicom图像
    
    filePath = strcat( outFolderPath,num2str(i*5+4),'.DCM');
    ImgInfoB.InstanceNumber = i * 5 + 5;
    ImgInfoB.SliceThickness = 1;
    ImgInfoB.Filename = filePath;
    ImgInfoB.SliceLocation = ImgInfoB.SliceLocation - 1;
    ImgInfoB.ImagePositionPatient(3,1) = ImgInfoB.ImagePositionPatient(3,1)-1;
    dicomwrite(Img4,filePath,ImgInfoB);%写入Dicom图像
    
    filePath = strcat( outFolderPath,num2str((i+1)*5+1),'.DCM');
    ImgInfoT.InstanceNumber = (startInstanceNumber + i) * 5 + 1;
    ImgInfoT.SliceThickness = 1;
    ImgInfoT.Filename = filePath;
    dicomwrite(ImgT,filePath,ImgInfoT);%写入Dicom图像 
end
CT图像三维重建是指通过医学诊断仪获取连续的二维切面图像,然后将这 些二维图像之间的位置和灰度信息输入计算机,在计算机上进行相应的组合和处理, 最后在显示器上再现人体该器官的立体影像并描绘出器官的三维图像信息。 CT图像三维重建能为医生显示具有真实感的三维图形,便于他们从多角度、 多次进行观察和分析。本文具体研究了医学图像三维重建所涉及的关键技术之一— —插值技术,提出了适合于 CT图像重建的轮廓形状插值算法,同时通过将该算 法应用到 CT图像三维重建系统中验证其合理性与有效性。 断图像插值是三维重建过程中的一个必要环节。本文对传统的插值方法进行了 归类,分为灰度插值、形状插值和小波插值三种,通过对这些方法中常用算法的分析 和比较,在实验的基础上讨论了各种算法的优缺点及适用范围,提出了一种基于轮廓 形状的 CT图像插值算法。 本文将 CT图像的三维重建系统划分为数据获取模块、数据预处理模块、图 像插值模块、绘制模块和显示模块,并使用相关技术实现了各模块的功能。在详细叙 述插值模块的设计与实现的基础上,将基于轮廓形状的插值算法和其它一些插值算法 应用于系统中进行比较。实验结果表明基于轮廓形状的插值算法能够比较好地保持图 像的边缘轮廓,有效地解决了梯田效应问题,为医学教学以及临床诊断治疗提供了比 较好的辅助。 论文最后对所作的工作进行了总结,并展望了下一步的研究工作。 关键词:CT图像,三维重建,插值技术,轮廓形状
if __name__ == '__main__': # -------------Adjustable global parameters---------- n=512 # pixel number m=10 # number of time phases angle = 5 # #sample points = 360/angle on the boundary numOfAngles = int(180/angle) numOfContourPts = int(360/angle) labelID = 1 # 勾画的RS文件中第几个轮廓为GTV # path of the input data folder = 'E:\\MedData\\4DCT-202305\\' #patient = '0007921948' # 缺少时间信息 patient = '0000726380' # 病人的编号 # 呼吸曲线数据文件 vxpPath = folder+patient+'\\0000726380\\0000726380_20230420_143723.vxp' # Save the generated figures to the latex file path figPath = "D:\\HUNNU\\Research\\DMD\\4D-CT\\latex-DMD插值\\modify202305\\figure\\" # -------------Auto generated global parameters---------- # 每个dicom文件包含多少横截面 name = os.listdir(folder+patient+'\\0') cuts = [] for i in range(len(name)): if 'CT' in name[i][0:2]: cuts.append(i+1) cuts = np.array(cuts) # phase name times = np.linspace(0,90,10) # image pixel coordinate nums = np.linspace(0,n-1,n) x,y = np.meshgrid(nums,nums) # 输出dicom头文件信息 filename = folder+patient+'\\0\\CT.{}'.format(patient)+'.Image 1.dcm' print('CT dicom file information:') info = loadFileInformation(filename) # 像素之间的间距,包括列间距和行间距,单位mm SliceThickness = info['SliceThickness'] # Z轴的扫描分辨率,单位mm pixelSpace = info['pixelSpace'] # 一个像素所占的实际体积 pixelVol = float(pixelSpace[0])*float(pixelSpace[0])*float(SliceThickness) print('sliceThickness=',SliceThickness,' pixelSpace=',pixelSpace)
06-02
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值