dicom文件的处理——pydicom库的使用

1. 计算两个切片的质量（dicom格式）

import cv2
from pydicom import dcmread
import numpy as np
from skimage.metrics import structural_similarity as compare_ssim
from skimage.metrics import peak_signal_noise_ratio as compare_psnr
from skimage.metrics import mean_squared_error as compare_mse

# 读取第一个DICOM文件
d1 = dcmread('PA1/SE0/IM0.IMA')
image_data1 = d1.pixel_array

# 读取第二个DICOM文件
d2 = dcmread('PA1/SE1/IM0.IMA')
image_data2 = d2.pixel_array

# 将图像转换为8位无符号整型格式
image_data1 = image_data1.astype(np.uint16)
image_data2 = image_data2.astype(np.uint16)

# 计算SSIM值
ssim_value = compare_ssim(image_data1, image_data2)
psnr_value = compare_psnr(image_data1,image_data2)
mse_value = compare_mse(image_data1,image_data2)
print(ssim_value, psnr_value, mse_value)

2. 对dicom文件进行重采样

之前的size大小不一，因此需要重采样到同一个size

import pydicom
import numpy as np

import pydicom
from PIL import Image
import os
from matplotlib import pyplot as plt
import lpips
from matplotlib import pyplot


def get_filelist(folder_path):
    file_names = os.listdir(folder_path)
    return file_names


def resample_dicom(ds, out_file, PA, idx):
    image_data = ds.pixel_array

    # 设置目标分辨率（这里以2x2作为示例）
    target_resolution = (512, 512)
    # 创建PIL图像对象
    image = Image.fromarray(image_data).resize(target_resolution)

    # 转换为numpy数组
    resampled_data = np.asanyarray(image)
    ds.PixelRepresentation = 1

    if ds.BitsAllocated == 16:
        resampled_data = resampled_data.astype(np.uint16, casting="safe")
    elif ds.BitsAllocated == 8:
        resampled_data = resampled_data.astype(np.int8, casting="safe")
    else:
        raise Exception("unknow Bits Allocated value in dicom header")

    # 设置新的像素数组
    ds.PixelData = resampled_data.tobytes()
    ds.Rows, ds.Columns = resampled_data.shape
    # 保存为新的DICOM文件
    # ds.save_as("output_file.dcm")
    ds.save_as("CTA-GAN1/"+ "PA"+str(PA+1)+"/"+out_file+"/IM"+str(idx)+".IMA")

for PA in range(10):

    SE0_path = "CTA-GAN/PA" + str(PA + 1) + "/SE0"
    SE1_path = "CTA-GAN/PA" + str(PA + 1) + "/SE1"

    files_name = get_filelist(SE0_path)
    for idx in range(len(files_name)):
        ds1 = pydicom.dcmread(SE0_path+"/IM"+str(idx)+".IMA",  force=True)  # 读取头文件
        ds2 = pydicom.dcmread(SE1_path+"/IM"+str(idx)+".IMA",  force=True)  # 读取头文件

        resample_dicom(ds1,"SE0", PA ,idx)
        resample_dicom(ds2, "SE1", PA, idx)

3. 创建训练、测试的txt，用于加载数据集

import os


os.remove('train.txt')
os.remove('val.txt')
os.remove('test.txt')

f1 = open("train.txt", "w")  # 564
f2 = open("val.txt", "w")  # 187
f3 = open("test.txt", "w")  # 188


def get_filelist(folder_path):
    file_names = os.listdir(folder_path)
    return file_names


def write_txt(path, txt):
    files_name = get_filelist(path)

    for i in range(len(files_name)):
        # aa = "data/" + path + "/" + files_name[i]
        aa = "data/" + path + "/" + "IM" + str(i) + ".IMA"
        txt.writelines(aa + "\n")


for i in range(8):
    SE0_path = "CTA-GAN/PA" + str(i + 1) + "/SE0"
    SE1_path = "CTA-GAN/PA" + str(i + 1) + "/SE1"

    write_txt(SE0_path, f1)
    write_txt(SE1_path, f1)
    

SE0_path = "CTA-GAN/PA9/SE0"
SE1_path = "CTA-GAN/PA9/SE1"

write_txt(SE0_path, f2)
write_txt(SE1_path, f2)

SE0_path = "CTA-GAN/PA10/SE0"
SE1_path = "CTA-GAN/PA10/SE1"

write_txt(SE0_path, f3)
write_txt(SE1_path, f3)

f1.close()
f1.close()
f2.close()

4. 对dicom进行rename

import os

def count_files(folder):
    file_count = len([name for name in os.listdir(folder) if os.path.isfile(os.path.join(folder, name))])
    return file_count

def get_filelist(folder_path):
    file_names = []
    for file in os.listdir(folder_path):
        if not os.path.isdir(os.path.join(folder_path, file)):
            file_names.append(file)
    return file_names


def rename_file(old_file, idx, file_path):
    # 定义原始文件路径和新文件名
    # old_file = "path/to/original.txt"
    old_file = os.path.join(file_path, old_file)
    new_name = os.path.join(file_path, "IM"+str(idx)+".IMA")

    try:
        # 调用rename()函数来重命名文件
        os.rename(old_file, new_name)
        print("文件已成功重命名为", new_name)
    except FileNotFoundError:
        print("未找到指定的文件")
    except Exception as e:
        print("发生错误：", str(e))



for i in range(10):
    SE0_path = "PA"+str(i+1)+"/SE0"
    SE1_path = "PA"+str(i+1)+"/SE1"
    num_file = count_files(SE0_path)

    files_name = get_filelist(SE0_path)
    for i in range(num_file):
        rename_file(files_name[i], i, SE0_path)

    files_name = get_filelist(SE1_path)
    for i in range(num_file):
        rename_file(files_name[i], i, SE1_path)

5. 删除文件夹下面的所有文件(只删除文件,不删除文件夹)

#删除文件夹下面的所有文件(只删除文件,不删除文件夹)
import os
import shutil

def del_file(path_data):
    for i in os.listdir(path_data) :# os.listdir(path_data)#返回一个列表，里面是当前目录下面的所有东西的相对路径
        file_data = path_data + "\\" + i#当前文件夹的下面的所有东西的绝对路径
        if os.path.isfile(file_data) == True:#os.path.isfile判断是否为文件,如果是文件,就删除.如果是文件夹.递归给del_file.
            os.remove(file_data)
        else:
            del_file(file_data)
path_data = r"CTA-GAN1"
del_file(path_data)

6. NLM滤波

import cv2
import numpy as np
import os

"""
基于相似性原理的图像去噪算法，与传统的局部均值滤波（平均，中值）不同，NLM 可以处理复杂的噪声模型，更好的保留图像细节
总得来说，NLM 算法的优缺点如下

优点：
去除噪音的同时，能够有效的保留图像的细节和纹理，降噪效果较好
不需要先验知识，可以处理多种不同类型的图像噪声
原理简单，易于实现

缺点：
复杂度高，运算速度较慢
对噪声强度和分布敏感，需要对参数进行合理设置才有好效果
处理低对比度和均匀纹理时，效果不佳

局部滤波（Local Filter）是使用当前像素的临近像素值来对图像进行降噪
非局部滤波（Non-local Filter）在这里的意思是寻找了相似块来进行降噪"""
def psnr(A, B):
    return 10*np.log(255*255.0/(((A.astype(float)-B)**2).mean()))/np.log(10)

def double2uint8(I, ratio=1.0):
    return np.clip(np.round(I*ratio), 0, 255).astype(np.uint8)

def make_kernel(f):
    """
    生成高斯核，越靠近中心位置的像素，权重越高
    """
    kernel = np.zeros((2*f+1, 2*f+1))
    for d in range(1, f+1):
        kernel[f-d:f+d+1, f-d:f+d+1] += (1.0/((2*d+1)**2))
    return kernel/kernel.sum()

def NLmeansfilter(I, h_=10, templateWindowSize=5,  searchWindowSize=11):
    f = templateWindowSize//2
    t = searchWindowSize//2
    height, width = I.shape[:2]
    padLength = t+f
    I2 = np.pad(I, padLength, 'symmetric')
    kernel = make_kernel(f)
    h = (h_**2)
    I_ = I2[padLength-f:padLength+f+height, padLength-f:padLength+f+width]

    average = np.zeros(I.shape)
    sweight = np.zeros(I.shape)
    wmax =  np.zeros(I.shape)
    for i in range(-t, t+1):
        for j in range(-t, t+1):
            if i==0 and j==0:
                continue
            I2_ = I2[padLength+i-f:padLength+i+f+height, padLength+j-f:padLength+j+f+width]
            # TODO: filter2D 的作用
            w = np.exp(-cv2.filter2D((I2_ - I_)**2, -1, kernel)/h)[f:f+height, f:f+width]
            sweight += w
            wmax = np.maximum(wmax, w)
            average += (w*I2_[f:f+height, f:f+width])
    # 原图像需要乘于最大权重参与计算
    average += (wmax*I)
    # sweight 为 weight 之和，用于计算 weight 的归一化
    sweight += wmax
    return average / sweight


def save_NLM(data_path):
    I = cv2.imread(data_path, 0)
    sigma = 20.0
    I1 = double2uint8(I + np.random.randn(*I.shape) * sigma)
    # cv2.imwrite('output_image0.jpg', I1)
    # print('噪声图像PSNR',psnr(I, I1))
    # R1  = cv2.medianBlur(I1, 5)
    # cv2.imwrite('output_image1.jpg', R1)
    # print('中值滤波PSNR',psnr(I, R1))
    # R2 = cv2.fastNlMeansDenoising(I1, None, sigma, 5, 11)
    # cv2.imwrite('output_image2.jpg', R2)
    # #print('opencv的NLM算法',psnr(I, R2))

    R3 = double2uint8(NLmeansfilter(I1.astype(float), sigma, 5, 11))
    cv2.imwrite(data_path, R3)
    # print('NLM PSNR', psnr(I, R3))




if __name__ == '__main__':
    """PSNR的数值越高,表示重建的图像质量越好。 PSNR的取值范围通常是0到∞,并且数值越大表示图像质量越好。
    具体来说,当PSNR接近∞时,意味着图像的质量损失非常小,几乎不可察觉;而当PSNR接近于0时,表示图像的质量损失非常大,重建的图像质量非常差。"""
    data_num = len(os.listdir("crop/png_ct"))


    for i in range(data_num):
        data_path = 'crop/png_ct/'+str(i)+'.png'
        save_NLM(data_path)

        data_path = 'crop/png_ct_aug/' + str(i+data_num) + '.png'
        save_NLM(data_path)

        # lable不用计算吧，因为就是二值的需要做任何处理
        """data_path = 'crop/png_label/' + str(i) + '.png'
        save_NLM(data_path)

        data_path = 'crop/png_label_aug/' + str(i+data_num) + '.png'
        save_NLM(data_path)"""

        """噪声图像PSNR 23.886751190851772
           中值滤波PSNR 27.526461361696413
           opencv的NLM算法 28.405511182856777
           NLM PSNR 28.618140873357483"""

7. 计算熵值

import cv2
from skimage.measure import shannon_entropy
import numpy as np

# 加载图像
img = cv2.imread('crop/png_ct/0.png', cv2.IMREAD_GRAYSCALE)

# 将图像转换为numpy数组
arr = np.array(img)

# 计算灰度直方图
histogram, _ = np.histogram(arr, bins=np.arange(0, 256), density=True)

# 计算熵值
entropy = -shannon_entropy(histogram)
print("Image entropy:", entropy)