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