LITS数据集预处理（二）

要求输出所有nii文件的size和spacing、origin、direction

size：nib库是（xyz),sitk为（zyx

origin：原始图像中心点在相机坐标系的位置

spacing：两个像素的间距

direction方向

写了一个代码可以同时输出文件名和对应spacing，size与其他可以同样代入

import os
import SimpleITK as itk
import numpy as np



def checksize(files):
    for i in (files):
        vol=itk.ReadImage(path+i)
        #filename,fileexp=os.path.split(i) 没用

        spacing=vol.GetSpacing()
        print(i,end=' ')   #执行不想换行操作,这里的i输出就是名字，如果需要去掉拓展名则需要用split操作
        print(spacing)


path='D:/test/'
files=os.listdir(path)
checksize(files)

原先预处理的代码被师兄改了一下，因为归一化是错的，还加了一些插值方法，这里记录一下：还有一个问题就是处理完后的数据比原来的大了很多，后来发现是因为有些图像spacing很大，重采样之后会大很多，有些很小会小很多。（这里我们重采样的方法是将spacing全部集合起来取中位数，大约取的是0.7，0.7,1.0）

import numpy as np
import SimpleITK as itk
import os
import matplotlib.pyplot as plt


def nii(files, newspacing=[0.7, 0.7, 1.0]):  # 设定新spacing为0.7,0.7,1.0
    for file in files:
        img = itk.ReadImage(path + file)  # 读取nii
        filename, fileexp = os.path.splitext(file)  # get filename
        original_spacing = img.GetSpacing()  # get spacing
        original_size = img.GetSize()  # get size

        new_size = [
            int(np.round(original_size[0] * original_spacing[0] / newspacing[0])),
            int(np.round(original_size[1] * original_spacing[1] / newspacing[1])),
            int(np.round(original_size[2] * original_spacing[2] / newspacing[2])),
        ]  # 新size

        resample = itk.ResampleImageFilter()
        resample.SetOutputDirection(img.GetDirection())
        resample.SetOutputOrigin(img.GetOrigin())
        resample.SetSize(new_size)
        resample.SetOutputSpacing(newspacing)

        if "segmentation" in file:     #for seg and img have different ways
            resample.SetInterpolator(itk.sitkNearestNeighbor)   #for seg,最临近插值，int8,否则0和1标记会直接插出0.5
            resample.SetOutputPixelType(itk.sitkUInt8)
        else:
            resample.SetInterpolator(itk.sitkLinear)
            resample.SetOutputPixelType(itk.sitkFloat32)  #for img线性插值，float32,
        img = resample.Execute(img)
        itk.WriteImage(img, r"H:/LITS/segmentation_processing/" + filename + ".nii")  # save and new files


path = r"H:/LITS/segmentation/"  # volume文件所在路径
files = os.listdir(path)  # 获取文件
nii(files)

接下来做本地化处理和导入网络预训练

 

6/12记录：今天接到了师兄的新任务，因为层厚太大了，所以想把多余的层截掉，然后再进行重采样，另外归一化也要想个办法实现，先记下来，明天去实现一下。（另外需要各向异性处理一下，一些spacing z太大的用最临近算法，比2.1小的用线性插值，原因是读图像的时候seg和vol的显示不流畅）

记录：d1、d2

6/13 各向异性处理

import numpy as np
import SimpleITK as itk
import os



def nii(files, newspacing=[0.7, 0.7, 1.0]):  # 设定新spacing为0.7,0.7,1.0
    for file in files:
        img = itk.ReadImage(path + file)  # 读取nii,记得要path+i才能够批量读取
        filename, fileexp = os.path.splitext(file)  # get filename
        original_spacing = img.GetSpacing()  # get spacing
        original_size = img.GetSize()  # get size

        new_size = [
            int(np.round(original_size[0] * original_spacing[0] / newspacing[0])),
            int(np.round(original_size[1] * original_spacing[1] / newspacing[1])),
            int(np.round(original_size[2] * original_spacing[2] / newspacing[2])),
        ]  # 新size

        resample = itk.ResampleImageFilter()
        resample.SetOutputDirection(img.GetDirection())
        resample.SetOutputOrigin(img.GetOrigin())
        resample.SetSize(new_size)
        resample.SetOutputSpacing(newspacing)

            #for seg and img have different ways
        resample.SetInterpolator(itk.sitkNearestNeighbor)   #for seg,最临近插值，int8,否则0和1标记会直接插出0.5
        resample.SetOutputPixelType(itk.sitkUInt8)

        # if original_spacing[2]>0.7*3 :   #for volume,>2.1最临近，<2.1线性
        #     resample.SetInterpolator(itk.sitkNearestNeighbor)
        #     resample.SetOutputPixelType(itk.sitkFloat32)  #for img线性插值，float32,
        # else:
        #     resample.SetInterpolator(itk.sitkLinear)
        #     resample.SetOutputPixelType(itk.sitkFloat32)
        img = resample.Execute(img)
        itk.WriteImage(img, "H:/MedicinePublicData/LITS2017/seg_pro/" + filename + ".nii")  # save and new files


path = "H:/MedicinePublicData/LITS2017/LITS2017/segmentation/"  # volume文件所在路径
files = os.listdir(path)  # 获取文件
nii(files)

修改了一下想法，觉得还是先重采样再crop比较好。crop是根据mask的位置上下裁的，这样数据比较小，训练方便。

之前没有做过，所以拿两个数据读入一下，看看能不能读到z轴的上下限

import os
import SimpleITK as itk
import numpy as np


def getsize(segfiles):
    for i in (segfiles):
        seg = itk.ReadImage(path2 + i)  # read seg
        seg_size = seg.GetSize()  # get seg size
        print(seg_size)

        seg_a=itk.GetArrayFromImage(seg)
        mask_index = np.argwhere(seg_a!= 0)  # return mask!=0 whole coords
        (d1,h1,w1),(d2,h2,w2)=np.max(mask_index,0),np.min(mask_index,0)
        d_max=d1
        d_min=d2
        print(d_max)
        print(d_min)

# def crop(volfiles,segfiles):
#     for i,j in (volfiles,segfiles):
#         vol=itk.ReadImage(path1+i) #read vol
#
#         filename1,fileexp1=os.path.split(i)  #get filename
#         filename2,fileexp2=os.path.split(j)
#         vol_size=vol.GetSize()  #get vol size
#         seg_size=seg.GetSize()  #get seg size
#
#         print(vol_size)
#         print(seg_size)
#
#         # mask_index=np.argwhere(seg!=0)  #return mask!=0 whole coords
#         # (d1,h1,w1),(d2,h2,w2)=np.max(mask_index,0),np.min(mask_index,0)

path1='H:/test/vol/'
path2='H:/test/label/'  #get segmentation
volfiles=os.listdir(path1)
segfiles=os.listdir(path2)
# crop(volfiles,segfiles)
getsize(segfiles)

可以读入，下一步是打算分别读图像和seg裁剪

存档今天的代码，报错了，明天继续改

import os
import SimpleITK as itk
import numpy as np


def getsize(segfiles,volfiles):
    for i,j in zip(segfiles,volfiles):
        seg = itk.ReadImage(path2 + i)  # read seg
        vol=itk.ReadImage(path1+j)
        filename1, fileexp1 = os.path.split(i)  # get filename
        filename2,fileexp2=os.path.split(j)
        seg_size = seg.GetSize()  # get seg size
        vol_size=vol.GetSize()
        print(seg_size)
        print(vol_size)


        seg_a=itk.GetArrayFromImage(seg)
        mask_index = np.argwhere(seg_a!= 0)  # return mask!=0 whole coords
        (d1,h1,w1),(d2,h2,w2)=np.max(mask_index,0),np.min(mask_index,0)
        d_max=d1 +3
        d_min=d2 -3

        newvol=vol[d_min:d_max,:,:]
        newseg=seg[d_min:d_max,:,:]

        itk.WriteImage(vol,'H:/test/voltest/'+filename1+'.nii')


# def crop(volfiles,segfiles):
#     for i,j in (volfiles,segfiles):
#         vol=itk.ReadImage(path1+i) #read vol
#
#         filename1,fileexp1=os.path.split(i)  #get filename
#         filename2,fileexp2=os.path.split(j)
#         vol_size=vol.GetSize()  #get vol size
#         seg_size=seg.GetSize()  #get seg size
#
#         print(vol_size)
#         print(seg_size)
#
#         # mask_index=np.argwhere(seg!=0)  #return mask!=0 whole coords
#         # (d1,h1,w1),(d2,h2,w2)=np.max(mask_index,0),np.min(mask_index,0)

path1='D:/test1/'
path2='H:/test/label/'  #get segmentation
volfiles=os.listdir(path1)
segfiles=os.listdir(path2)
# crop(volfiles,segfiles)
getsize(segfiles,volfiles)

6/17解决了报错的问题：取名字的函数写错了，晕

源代码不见了，拿隔壁的数据集的crop代码来，大概lits是要替换成z轴即可，crop也crop z 轴

import os
import SimpleITK as itk
import numpy as np


def crop(segfiles,volfiles):
    for i,j in zip(segfiles,volfiles):
        seg = itk.ReadImage(path2 + i)  # read seg
        vol=itk.ReadImage(path1+j)     #read vol
        filename1, fileexp1 = os.path.splitext(i)  # get filename:filename1:seg
        filename2,fileexp2=os.path.splitext(j)  #filename2:vol
        seg_size = seg.GetSize()  # get seg size
        vol_size=vol.GetSize()   #get vol size
        # print(seg_size)
        # print(vol_size)


        seg_a=itk.GetArrayFromImage(seg)
        mask_index = np.argwhere(seg_a!= 0)  # return mask!=0 whole coords
        (d1,h1,w1),(d2,h2,w2)=np.max(mask_index,0),np.min(mask_index,0)  #zyx
        h_max=h1 +3   #side
        h_min=h2 -3

        newvol=vol[:,h_min:h_max,:]  #new vol,对于3D图像可以这样直接CROP xyz
        newseg=seg[:,h_min:h_max,:]  #new seg

        itk.WriteImage(newvol,'H:/MedicinePublicData/Task02_Heart/imgtrain_cropy/'+filename2+'.nii')  #save vol
        itk.WriteImage(newseg, 'H:/MedicinePublicData/Task02_Heart/label_cropy/' + filename1 + '.nii')  #save seg


# def crop(volfiles,segfiles):
#     for i,j in (volfiles,segfiles):
#         vol=itk.ReadImage(path1+i) #read vol
#
#         filename1,fileexp1=os.path.split(i)  #get filename
#         filename2,fileexp2=os.path.split(j)
#         vol_size=vol.GetSize()  #get vol size
#         seg_size=seg.GetSize()  #get seg size
#
#         print(vol_size)
#         print(seg_size)
#
#         # mask_index=np.argwhere(seg!=0)  #return mask!=0 whole coords
#         # (d1,h1,w1),(d2,h2,w2)=np.max(mask_index,0),np.min(mask_index,0)

path1='H:/MedicinePublicData/Task02_Heart/imgtrain_pro/'   #get vol
path2='H:/MedicinePublicData/Task02_Heart/label_pro/'  #get segmentation
volfiles=os.listdir(path1)
segfiles=os.listdir(path2)
# crop(volfiles,segfiles)
crop(segfiles,volfiles)