itk和sitk的一些应用整理

1.sitk读取dicom并打印dicom中的tag

# 通过sitk读取dicom文件
def read_dicom_file_with_sitk(file_path, print_info=0):
    dcmImage = sitk.ReadImage(file_path, sitk.sitkUInt16)
    metadata = dcmImage.GetMetaDataKeys()
    for key in metadata:
        value = dcmImage.GetMetaData(key)
        print(f"{key}: {value}")
    if print_info == 1:
        print(type(dcmImage))
        print(dcmImage.GetSize())
        print(dcmImage.GetOrigin())
        print(dcmImage.GetSpacing())
        print(dcmImage.GetDirection())
        # print(dcmImage.GetMetaData('0002|0010')) # TransferSyntaxUID
        print(dcmImage.GetMetaData('0008|0060'))  # Modality
        print(dcmImage.GetMetaData('0010|0010'))  # PatientName
        print(dcmImage.GetMetaData('0028|0002'))  # SamplesPerPixel
        print(dcmImage.GetMetaData('0028|0100'))  # BitsAllocated
    window_center = float(dcmImage.GetMetaData('0028|1050'))  # WindowCenter
    window_width = float(dcmImage.GetMetaData('0028|1051'))  # WindowWidth

    image_array = sitk.GetArrayFromImage(dcmImage)[0]
    print(type(image_array), np.max(image_array), np.min(image_array))
    plt.imshow(image_array, cmap='gray')
    plt.axis('off')
    plt.show()
    return image_array, window_center, window_width

2.itk读取dicom文件夹并打印dicom 的tag信息

# 打印dicom的tag信息
def print_tag_info(metadata):
    # Print the key value pairs from the metadadictionary
    tag_keys = metadata.GetKeys()
    for tag_key in tag_keys:
        # Note the [] operator for the key
        try:
            label = itk.GDCMImageIO.GetLabelFromTag(tag_key, "")
            tagvalue = metadata[tag_key]
            print(label[1] + '(' + tag_key + ")=" + str(tagvalue))
        except RuntimeError:
            # Cannot pass specialized values into metadata dictionary.
            print("Cannot pass specialized value " + tag_key + " into metadadictionary")

    # 如何使用给定的标签来获取标签
    entryID = "0028|0102"
    if not metadata.HasKey(entryID):
        print("tag: " + entryID + " not found in series")
    else:
        # The second parameter is mandatory in python to get the
        # string label value
        label = itk.GDCMImageIO.GetLabelFromTag(entryID, "")
        tag_value = metadata[entryID]
        print(label[1] + " (" + entryID + ") is: " + str(tag_value))


# 通过itk读取dicom文件夹
def read_dicom_dir(dicom_dir, save_data=0, print_tag=0):
    PixelType = itk.ctype("signed short")
    Dimension = 3
    ImageType = itk.Image[PixelType, Dimension]

    namesGenerator = itk.GDCMSeriesFileNames.New()
    namesGenerator.SetUseSeriesDetails(True)
    namesGenerator.AddSeriesRestriction("0008|0021")
    namesGenerator.SetGlobalWarningDisplay(False)
    namesGenerator.SetDirectory(dicom_dir)

    seriesUID = namesGenerator.GetSeriesUIDs()

    if len(seriesUID) < 1:
        print("No DICOMs in: " + dicom_dir)
        sys.exit(1)

    print("The directory: " + dicom_dir)
    print("Contains the following DICOM Series: ")
    for uid in seriesUID:
        print(uid)

    for uid in seriesUID:
        seriesIdentifier = uid
        print("Reading: " + seriesIdentifier)
        fileNames = namesGenerator.GetFileNames(seriesIdentifier)

        reader = itk.ImageSeriesReader[ImageType].New()
        dicomIO = itk.GDCMImageIO.New()
        reader.SetImageIO(dicomIO)
        reader.SetFileNames(fileNames)
        reader.ForceOrthogonalDirectionOff()
        reader.Update()
        if save_data == 1:
            writer = itk.ImageFileWriter[ImageType].New()
            outFileName = os.path.join(dicom_dir, seriesIdentifier + ".nrrd")
            writer.SetFileName(outFileName)
            writer.UseCompressionOn()
            writer.SetInput(reader.GetOutput())
            print("Writing: " + outFileName)
            writer.Update()
        outItkImage = reader.GetOutput()
        print('origin:', outItkImage.GetOrigin())
        vtkImage = itk.vtk_image_from_image(outItkImage)
        dims = vtkImage.GetDimensions()
        spacings = vtkImage.GetSpacing()
        origin = vtkImage.GetOrigin()
        center = vtkImage.GetCenter()
        scalarRange = vtkImage.GetScalarRange()
        print(dims, spacings, origin, center)
        print(scalarRange)
        metadata = dicomIO.GetMetaDataDictionary()
        if print_tag:
            print_tag_info(metadata)
        return vtkImage, outItkImage

3.itk读取单个nii.gz文件和多个nii.gz文件合并

# 读取nii.gz
def read_niigz_file(file_path):
    PixelType = itk.SS
    Dimension = 3
    ImageType = itk.Image[PixelType, Dimension]
    reader = itk.ImageFileReader[ImageType].New()
    reader.SetFileName(file_path)
    reader.Update()
    image = reader.GetOutput()
    return itk.vtk_image_from_image(image), image


# 读取多个nii.gz文件,并将其合并为一个nii.gz
def read_niigz_files(file_paths):
    sum_itk_image = None
    for file_path in file_paths:
        vtk_image, itk_image = read_niigz_file(file_path)
        image_array = itk.GetArrayFromImage(itk_image)
        if sum_itk_image is None:
            sum_itk_image = image_array
        else:
            sum_itk_image = sum_itk_image + image_array

    res_itk_image = itk.GetImageFromArray(sum_itk_image)
    res_vtk_image = itk.vtk_image_from_image(res_itk_image)
    return res_vtk_image, res_itk_image

4.itk读取nrrd文件 (翻转itkImage)和多个nrrd文件合并

# 翻转一个itkImage
def flip_itk_image(in_image, flip_axes=[False, False, False]):
    flip_image = itk.FlipImageFilter[type(in_image)].New()
    flip_image.SetFlipAxes(flip_axes)
    flip_image.SetInput(in_image)
    flip_image.Update()


# 通过itk读取nrrd文件
def read_nrrd(file_name):
    PixelType = itk.ctype("unsigned char")
    Dimension = 3
    ImageType = itk.Image[PixelType, Dimension]
    reader = itk.ImageFileReader[ImageType].New()
    reader.SetFileName(file_name)
    reader.Update()

    out_itk_image = reader.GetOutput()
    image_origin = out_itk_image.GetOrigin()
    print('out_itk_image origin:', image_origin)
    vtkImage = itk.vtk_image_from_image(out_itk_image)
    vtkImage.SetOrigin(image_origin)
    print('spacing:', vtkImage.GetSpacing())
    print('dims:', vtkImage.GetDimensions())
    return vtkImage, out_itk_image


# 读取多个nrrd文件并合并为一个文件
def read_nrrd_files(file_paths):
    sum_itk_image = None
    for file_path in file_paths:
        vtk_image, itk_image = read_nrrd(file_path)
        image_array = itk.GetArrayFromImage(itk_image)
        if sum_itk_image is None:
            sum_itk_image = image_array
        else:
            sum_itk_image = sum_itk_image + image_array

    res_itk_image = itk.GetImageFromArray(sum_itk_image)
    res_vtk_image = itk.vtk_image_from_image(res_itk_image)
    return res_vtk_image, res_itk_image



def sum_nrrd_file(base_dir):
    all_files = []
    for root, dirs, files in os.walk(base_dir):
        for file in files:
            if file[-4:] == 'nrrd':
                all_files.append(os.path.join(root, file))

    data1 = read_nrrd(all_files[0])
    for file in all_files:
        tmp_data = read_nrrd(file)
        indices = np.where(tmp_data > 0)
        data1[indices] = tmp_data[indices]

    print('read end...')

    image = sitk.GetImageFromArray(data1)
    image.SetSpacing([0.732422, 0.732422, 1.0])
    image.SetOrigin([-187.134, -365.134, -1727])
    sitk.WriteImage(image, "all.nrrd")
    print('end...')

5.更改itkImage的类型(获取itkImage的最大最小值)

# 获取itkImage中的最大值和最小值
def get_image_min_max_value(in_image):
    ImageCalculatorFilterType = itk.MinimumMaximumImageCalculator[type(in_image)]
    imageCalculatorFilter = ImageCalculatorFilterType.New()
    imageCalculatorFilter.SetImage(in_image)
    imageCalculatorFilter.Compute()
    minimumLocation = imageCalculatorFilter.GetIndexOfMinimum()
    maximumLocation = imageCalculatorFilter.GetIndexOfMaximum()
    minimumValue = imageCalculatorFilter.GetMinimum()
    maximumValue = imageCalculatorFilter.GetMaximum()
    print('minimumValue:', minimumValue)
    print('maximumValue:', maximumValue)
    return minimumValue, maximumValue


# 打印itkImage中的一些信息
def print_data_info(itk_image):
    imOrigin = itk_image.GetOrigin()
    imRes = itk_image.GetSpacing()
    imRegion = itk_image.GetBufferedRegion()
    imSize = imRegion.GetSize()
    print(imSize, imRes, imOrigin)


# 转换image的类型为unsigned short
def change_image_type(in_image):
    min_value, max_value = get_image_min_max_value(in_image)
    rescaler = itk.RescaleIntensityImageFilter[type(in_image), itk.Image[itk.US, 3]].New()
    rescaler.SetInput(in_image)
    rescaler.SetOutputMinimum(0)
    rescaler.SetOutputMaximum(max_value-min_value)
    rescaler.Update()
    return rescaler.GetOutput()

6.保存nrrd文件为mhd文件

# 保存vtkImageData为mhd文件
def save_imagedata_mhd(imagedata, file_path, file_path_raw):
    writer = vtk.vtkMetaImageWriter()
    writer.SetInputData(imagedata)
    writer.SetFileName(file_path)
    writer.SetRAWFileName(file_path_raw)
    writer.Write()
    print('save mhd succeed...')


# 获取体的一部分
def get_part_imagedata(imagedata, voi):
    extractVOI = vtk.vtkExtractVOI()
    extractVOI.SetInputData(imagedata)
    extractVOI.SetVOI(voi[0], voi[1], voi[2], voi[3], voi[4], voi[5])
    extractVOI.Update()
    save_file_name = 'part_voi.mhd'
    save_file_name_raw = 'part_voi.raw'
    save_imagedata_mhd(extractVOI.GetOutput(), save_file_name, save_file_name_raw)