dicom影像坐标转换

木鱼敲僧

已于 2025-01-16 15:21:02 修改

阅读量329

点赞数 10

文章标签：计算机视觉人工智能 python

于 2025-01-16 15:16:15 首次发布

本文链接：https://blog.youkuaiyun.com/qq_44942100/article/details/145174051

版权

有个需求要把病人坐标转换到图像像素坐标上，然后描点、填充生成mask。

搞了两天头都快秃了，试了网上的好多方法都对不上，还是要结合实际情况多试试。

1、roi提取关键信息为csv

import os
import csv
import re

def extract_points_with_name_from_roi(roi_file_path, csv_file_path):
    data_points = []
    current_name = None
    name_pattern = re.compile(r'name:\s*(\S+)')

    with open(roi_file_path, 'r') as f:
        read_points = False
        for line in f:
            line = line.strip()

            if line.startswith('name:'):
                match = name_pattern.search(line)
                if match:
                    current_name = match.group(1)

            elif 'points={' in line:
                read_points = True
                continue

            elif '};' in line and read_points:
                read_points = False
                continue

            elif read_points:
                if line:
                    values = line.split()
                    if len(values) >= 3:
                        try:
                            x, y, z = float(values[0]), float(values[1]), float(values[2])
                            data_points.append([current_name, x, y, z])
                        except ValueError:
                            pass

    with open(csv_file_path, 'w', newline='') as csvfile:
        csv_writer = csv.writer(csvfile)
        csv_writer.writerow(['name', 'x', 'y', 'z'])
        csv_writer.writerows(data_points)

source_dir = './npc'

for patient_folder in os.listdir(source_dir):
    patient_folder_path = os.path.join(source_dir, patient_folder)
    if os.path.isdir(patient_folder_path):
        for file_name in os.listdir(patient_folder_path):
            if file_name.endswith('.roi'):
                roi_file_path = os.path.join(patient_folder_path, file_name)
                csv_file_name = file_name.replace('.roi', '.csv')
                csv_file_path = os.path.join(patient_folder_path, csv_file_name)
                extract_points_with_name_from_roi(roi_file_path, csv_file_path)

根据z唯一值数量（影像数量）重命名csv：

import os
import pandas as pd

source_dir = './npc'

for patient_folder in os.listdir(source_dir):
    patient_folder_path = os.path.join(source_dir, patient_folder)
    if os.path.isdir(patient_folder_path):
        for file_name in os.listdir(patient_folder_path):
            if file_name.endswith('.csv'):
                csv_file_path = os.path.join(patient_folder_path, file_name)
                df = pd.read_csv(csv_file_path)
                
                if 'z' in df.columns:
                    unique_z_count = df['z'].nunique()
                    new_file_name = file_name.replace('.csv', f'_{unique_z_count}.csv')
                    new_file_path = os.path.join(patient_folder_path, new_file_name)
                    
                    os.rename(csv_file_path, new_file_path)

删除不满足影像数量的csv：

import os
import pydicom

def count_dicom_files(dicom_dir):
    return len([file for file in os.listdir(dicom_dir) if file.endswith('.dcm')])

def process_patient_folders(base_dir):
    for folder in os.listdir(base_dir):
        folder_path = os.path.join(base_dir, folder)
        if not os.path.isdir(folder_path):
            continue

        dicom_dir = os.path.join(folder_path, 'ImageSet_0.DICOM')
        if not os.path.exists(dicom_dir):
            print(f"DICOM directory not found for {folder}. Skipping...")
            continue

        dicom_count = count_dicom_files(dicom_dir)
        for file in os.listdir(folder_path):
            if file.endswith('.csv'):
                csv_path = os.path.join(folder_path, file)
                try:
                    count_in_filename = int(file.split('_')[-1].split('.')[0])
                    if count_in_filename != dicom_count:
                        os.remove(csv_path)
                except ValueError:
                    print(f"Invalid format in file name: {file}. Skipping...")

base_dir = './npc'
process_patient_folders(base_dir)

2、坐标转换

首先观察z值，发现唯一值正好就是图像的数量，而且和图像的SliceLocation是对应的，有个统一的倍数关系。

import numpy as np

data = np.genfromtxt('./demo.csv', delimiter=',', skip_header=1)
unique_depth_values = np.unique(data[:, 3])
num_unique_values = len(unique_depth_values)

print(num_unique_values)
for value in unique_depth_values:
    print(value)

然后就是x和y了，确定单位（厘米or 毫米），x和y分别除以像素间距（单位毫米）pixel_spacing[0]和pixel_spacing[1]，就能得到偏移坐标。

最关键的是确定原点的位置，目前还没找到准确的找原点方式。

读取所有影像尺寸的唯一值：

import os
import pydicom

def get_dicom_image_size(file_path):
    try:
        ds = pydicom.dcmread(file_path)
        return ds.Rows, ds.Columns
    except Exception as e:
        print(f"Error reading {file_path}: {e}")
        return None

def traverse_and_collect_unique_sizes(root_dir):
    unique_sizes = set()
    for subdir, _, files in os.walk(root_dir):
        for file in files:
            if file.lower().endswith('.dcm'):
                file_path = os.path.join(subdir, file)
                size = get_dicom_image_size(file_path)
                if size:
                    unique_sizes.add(size)
    return sorted(unique_sizes)

root_directory = "./npc"
unique_sizes = traverse_and_collect_unique_sizes(root_directory)

print("Number of unique image sizes: ", len(unique_sizes))
print("Unique image sizes: ", unique_sizes)

通过结果反推发现有一部分图像的原点是图像尺寸的正中心，比如512x512的图像，那么原点就是(256,256)，生成的结果感觉上是准的，后面有机会弄清楚了再修改吧。

import os
import pandas as pd
import pydicom

def get_first_dicom_file(dicom_dir):
    for file in os.listdir(dicom_dir):
        file_path = os.path.join(dicom_dir, file)
        if os.path.isfile(file_path) and file.lower().endswith('.dcm'):
            return file_path
    raise FileNotFoundError(f"No valid DICOM files found in {dicom_dir}")

def process_csv(csv_file, dicom_dir):
    df = pd.read_csv(csv_file)
    dicom_file = get_first_dicom_file(dicom_dir)
    ds = pydicom.dcmread(dicom_file)
    pixel_spacing = [ps / 10 for ps in ds.PixelSpacing]
    z_unique_sorted = {val: rank + 1 for rank, val in enumerate(sorted(df['z'].unique(), reverse=True))}
    df['z'] = df['z'].map(z_unique_sorted)
    
    df['x'] = (256 + (df['x'] / pixel_spacing[0])).round().astype(int)
    df['y'] = (256 - (df['y'] / pixel_spacing[1])).round().astype(int)
    df['z'] = df['z'].round().astype(int)
    
    new_csv_file = os.path.join(os.path.dirname(csv_file), f"T_{os.path.basename(csv_file)}")
    df.to_csv(new_csv_file, index=False)

base_dir = './T'

for patient_folder in os.listdir(base_dir):
    patient_path = os.path.join(base_dir, patient_folder)
    if os.path.isdir(patient_path):
        dicom_dir = os.path.join(patient_path, 'ImageSet_0.DICOM')
        if os.path.isdir(dicom_dir):
            for file in os.listdir(patient_path):
                if file.endswith('.csv'):
                    csv_file = os.path.join(patient_path, file)
                    process_csv(csv_file, dicom_dir)

还有一部分图像如果按照256算的话，y坐标是对不上的，尝试过后发现可能是第一张切片位置ImagePositionPatient属性的y，那么代码就变成了：

import os
import pandas as pd
import pydicom

def get_first_dicom_file(dicom_dir):
    for file in os.listdir(dicom_dir):
        file_path = os.path.join(dicom_dir, file)
        if os.path.isfile(file_path) and file.lower().endswith('.dcm'):
            return file_path
    raise FileNotFoundError(f"No valid DICOM files found in {dicom_dir}")

def process_csv(csv_file, dicom_dir):
    df = pd.read_csv(csv_file)
    dicom_file = get_first_dicom_file(dicom_dir)
    ds = pydicom.dcmread(dicom_file)
    
    pixel_spacing = [ps / 10 for ps in ds.PixelSpacing]
    image_position_y = abs(ds.ImagePositionPatient[1])
    
    z_unique_sorted = {val: rank + 1 for rank, val in enumerate(sorted(df['z'].unique(), reverse=True))}
    df['z'] = df['z'].map(z_unique_sorted)
    
    df['x'] = (256 + (df['x'] / pixel_spacing[0])).round().astype(int)
    df['y'] = (image_position_y - (df['y'] / pixel_spacing[1])).round().astype(int)
    df['z'] = df['z'].round().astype(int)
    
    new_csv_file = os.path.join(os.path.dirname(csv_file), f"T_{os.path.basename(csv_file)}")
    df.to_csv(new_csv_file, index=False)

base_dir = './F'

for patient_folder in os.listdir(base_dir):
    patient_path = os.path.join(base_dir, patient_folder)
    if os.path.isdir(patient_path):
        dicom_dir = os.path.join(patient_path, 'ImageSet_0.DICOM')
        if os.path.isdir(dicom_dir):
            for file in os.listdir(patient_path):
                if file.endswith('.csv'):
                    csv_file = os.path.join(patient_path, file)
                    process_csv(csv_file, dicom_dir)

到这里后我突然觉得读取ImagePositionPatient属性里的x,y才是准确的，可生成的点还是对不上，真是奇了怪了......此刻想要获取知识的心情达到了顶峰。

3、描点

import os
import numpy as np
import nibabel as nib
import csv
from collections import defaultdict

def generate_nii_from_csv(data, output_nii_path, image_size, depth):
    nii_array = np.zeros((image_size, image_size, depth), dtype=np.uint8)
    for r, c, d in data:
        d = d - 1  # 将切片索引从1开始调整为0基索引
        if 0 <= r < image_size and 0 <= c < image_size and 0 <= d < depth:
            nii_array[r, c, d] = 1
    nii_image = nib.Nifti1Image(nii_array, affine=np.eye(4))
    nib.save(nii_image, output_nii_path)

def process_patient_folders(base_dir, image_size=512):
    for folder in os.listdir(base_dir):
        folder_path = os.path.join(base_dir, folder)
        if not os.path.isdir(folder_path):
            continue
        
        data_dict = defaultdict(list)
        
        for file in os.listdir(folder_path):
            if file.startswith('T_') and file.endswith('.csv'):
                csv_file_path = os.path.join(folder_path, file)
                try:
                    depth = int(file.split('_')[-1].split('.')[0])
                except ValueError:
                    continue
                
                with open(csv_file_path, 'r') as csvfile:
                    csv_reader = csv.reader(csvfile)
                    next(csv_reader)
                    for row in csv_reader:
                        name, r, c, d = row
                        r, c, d = map(int, [r, c, d])
                        data_dict[name].append((r, c, d))
        
        mask_dir = os.path.join(folder_path, 'curve')
        if not os.path.exists(mask_dir):
            os.makedirs(mask_dir)
        
        for name, data in data_dict.items():
            output_nii_path = os.path.join(mask_dir, f"{name}.nii.gz")
            generate_nii_from_csv(data, output_nii_path, image_size, depth)

base_dir = './F'
process_patient_folders(base_dir)

4、填充

import os
import numpy as np
import nibabel as nib
import csv
from skimage import measure
from collections import defaultdict
from scipy.ndimage import binary_fill_holes
from skimage.draw import polygon

def generate_nii_from_csv(data, output_nii_path, image_size, depth):
    nii_array = np.zeros((image_size, image_size, depth), dtype=np.uint8)
    for d in range(depth):
        slice_points = [(r, c) for r, c, z in data if z == d + 1]  # 调整索引从1开始
        if len(slice_points) > 2:
            rows, cols = zip(*slice_points)
            rr, cc = polygon(rows, cols, (image_size, image_size))
            labels = np.zeros_like(nii_array[:, :, d])
            labels[rr, cc] = 1
            labeled_array = measure.label(labels)
            for region in measure.regionprops(labeled_array):
                if region.area > 2:
                    region_coords = region.coords
                    region_rr, region_cc = region_coords[:, 0], region_coords[:, 1]
                    nii_array[region_rr, region_cc, d] = 1
            nii_array[:, :, d] = binary_fill_holes(nii_array[:, :, d]).astype(np.uint8)
    
    nii_image = nib.Nifti1Image(nii_array, affine=np.eye(4))
    nib.save(nii_image, output_nii_path)

def process_patient_folders(base_dir, image_size=512):
    for folder in os.listdir(base_dir):
        folder_path = os.path.join(base_dir, folder)
        if not os.path.isdir(folder_path):
            continue
        
        data_dict = defaultdict(list)
        
        for file in os.listdir(folder_path):
            if file.startswith('T_') and file.endswith('.csv'):
                csv_file_path = os.path.join(folder_path, file)
                try:
                    depth = int(file.split('_')[-1].split('.')[0])
                except ValueError:
                    continue
                
                with open(csv_file_path, 'r') as csvfile:
                    csv_reader = csv.reader(csvfile)
                    next(csv_reader)
                    for row in csv_reader:
                        name, r, c, d = row
                        r, c, d = map(int, [r, c, d])
                        data_dict[name].append((r, c, d))
        
        mask_dir = os.path.join(folder_path, 'mask')
        if not os.path.exists(mask_dir):
            os.makedirs(mask_dir)
        
        for name, data in data_dict.items():
            output_nii_path = os.path.join(mask_dir, f"{name}.nii.gz")
            generate_nii_from_csv(data, output_nii_path, image_size, depth)

base_dir = './F'
process_patient_folders(base_dir)

附终端删除无用文件的命令：