首页 > 其他分享 >dicom影像坐标转换

dicom影像坐标转换

时间:2025-01-16 15:58:10浏览次数:3  
标签:dicom dir 坐标 file path csv os 影像

有个需求要把病人坐标转换到图像像素坐标上,然后描点、填充生成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)

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

find ./npc/nii -type d -name ImageSet_0.DICOM -exec rm -r {} +

find ./npc/nii -type f -name "*.csv" -exec rm -f {} +

标签:dicom,dir,坐标,file,path,csv,os,影像
From: https://blog.csdn.net/qq_44942100/article/details/145174051

相关文章

  • 【2024遥感应用组一等奖】基于ENVI Deep Learning的多源SAR影像海带水淹区域提取与分
    作品介绍01应用背景近年来,随着全球气候变化的加剧,海岸带地区面临的洪水威胁日益严重。这些极端天气事件,特别是由热带气旋引发的洪水,不仅给沿海地区带来了巨大的经济损失,还对当地居民的安全和生活产生了深远影响。因此,能够准确识别海岸地区的水淹区域对于防灾减灾、制定应对策......
  • 【深度学习地学应用|滑坡制图、变化检测、多目标域适应、感知学习、深度学习】跨域大
    【深度学习地学应用|滑坡制图、变化检测、多目标域适应、感知学习、深度学习】跨域大尺度遥感影像滑坡制图方法:基于原型引导的领域感知渐进表示学习(五)【深度学习地学应用|滑坡制图、变化检测、多目标域适应、感知学习、深度学习】跨域大尺度遥感影像滑坡制图方法:基于原型引......
  • 基于YOLOv5的医学影像病变区域识别:深度学习在医学诊断中的应用
    随着医疗技术的进步,医学影像成为了辅助医生进行疾病诊断的重要工具。医学影像不仅能帮助医生观察到患者体内的病变区域,还能为疾病的早期发现和精准治疗提供关键线索。传统的医学影像分析方法依赖于医生的经验和人工判断,效率低且容易受到人为因素的影响。而随着深度学习技术,......
  • python脚本实现经纬度和大地高与ECEF坐标互转
    importmath#地心地固坐标系(ECEF)转经纬度和大地高defecef2lla(x,y,z):#初始近似d=0for_inrange(32):#最大迭代次数设为32,可根据实际情况调整#计算临时变量R_prime=math.sqrt(x**2+y**2+(z-d)**2)......
  • ros2笔记-5.3 C++中地图坐标系变换
    本节继续跟小鱼老师学习5.3.需求背景:地图坐标系为map,机器人坐标系为baselink,目标点为target_point,已知map到baselink之间的关系,map到target_point关系。要控制机器人到达目标点,就需要知道目标点target_point和机器人base_link之间的关系。5.3.1通过C++发布静态TF目标......
  • 【ENVI初学】查看影像空间分辨率及重采样
    一、查看影像空间分辨率方法一:查看影像的元数据在ENVI中打开影像后右键选择ViewMetadata打开元数据窗口 MapInfo中的PixelSize指示了每个像素的地理跨度(单位通常是地理坐标系统中的度、投影坐标系统中的米、千米或其他投影单位)例如,当PixelsizeX=0.05 时1. 如果......
  • 前端学习openLayers配合vue3(两个坐标之间线的绘制)
    上节我们学了点的绘制,今天我们来学习一下线的绘制关键代码constlineLayer=newVectorLayer({source:newVectorSource(),});map.addLayer(lineLayer);letfeature=newFeature({//北京到上海的经纬度geometry:newLineString([[116.46,39.92],......
  • 请说说地图相关的坐标系有哪些?
    在地图相关领域中,坐标系是确定和定位地理要素(如点、线、面)位置的基础。对于前端开发而言,了解这些坐标系对于实现地图的可视化和交互至关重要。以下是一些常见的地图相关坐标系:地理坐标系:使用经度和纬度来表示地理位置。地球表面上的任意一点都可以通过一组经纬度坐标来唯一确......
  • 各坐标系是如何转换的?
    在前端开发中,坐标系转换是一个常见的需求,尤其是当处理地图数据或集成不同地图服务时。以下是一些常见的坐标系及其之间的转换方法:1.坐标系概述GCJ-02(火星坐标系):由中国国家测绘地理信息局制定,用于国内地图和位置服务,是一种加密坐标系。BD-09(百度坐标系):基于GCJ-02坐标系,由百度......
  • 赋能百强,霄云科技助力复旦大学附属肿瘤医院影像数据存储建设
    随着医疗技术的飞速发展,医学影像数据已成为医疗信息领域的核心组成部分,其数据量庞大且增长迅速,对存储、管理和使用提出了极高的要求。复旦大学附属肿瘤医院是我国成立最早、集医、教、研、防为一体的三级甲等肿瘤专科医院,年门诊量超200万人次、影像检查年数据增量超120TB。......