放射科的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