如何修改CT影像中心脏区域的Hounsfield Unit值为固定值

我拥有DICOM格式的CT影像与RT结构文件，希望将心脏区域的Hounsfield Unit（HU值）修改为固定值100。目前我已实现心脏HU值的计算与展示，但无法完成修改操作，恳请提供指导。

用户原实现代码：

import os
import pydicom
import numpy as np
from skimage.draw import polygon
import matplotlib.pyplot as plt

# Directory containing CT images
ct_directory = r"D:\ct directory"

# Path to your RT Structure file
rt_struct_path = r"D:\RTS2.dcm"

# Load CT images
ct_files = [f for f in os.listdir(ct_directory) if f.endswith('.dcm')]
ct_images = [pydicom.dcmread(os.path.join(ct_directory, f)) for f in ct_files]
ct_images.sort(key=lambda x: float(x.ImagePositionPatient[2]))

# Load RT Structure
rt_struct = pydicom.dcmread(rt_struct_path)

# Identify the specific ROI number for the heart (adjust 'heart' if necessary based on your DICOM  files)
heart_roi_number = None
for roi in rt_struct.StructureSetROISequence:
    if 'heart' in roi.ROIName.lower():  # Adjust the condition based on the exact name
        heart_roi_number = roi.ROINumber
        break

if heart_roi_number is None:
    raise ValueError("Heart ROI not found")

# Loop through each CT slice and process contours to calculate Hounsfield Units
for img in ct_images:
    try:
        img.contour_points = []  # Ensure this attribute is reset for each image

        # Check if this slice contains the heart contour
        heart_contour_found = False
        for roi_contour in rt_struct.ROIContourSequence:
            if roi_contour.ReferencedROINumber == heart_roi_number:
                for contour in roi_contour.ContourSequence:
                    contour_data = np.array(contour.ContourData).reshape((-1, 3))

                    # Process each point in the contour
                    for point in contour_data:
                        z_value = point[2]

                        # Find the corresponding CT slice for this point
                        if abs(img.ImagePositionPatient[2] - z_value) < 1e-2:

                            # Correctly convert this point's (x, y) to pixel coordinates
                            x, y = point[:2]
                            col = np.round((x - img.ImagePositionPatient[0]) /    img.PixelSpacing[0]).astype(int)
                            row = np.round((y - img.ImagePositionPatient[1]) /  img.PixelSpacing[1]).astype(int)

                            img.contour_points.append((row, col))

                            # If we find a slice with the heart contour, set heart_contour_found to True
                            if not heart_contour_found:
                                heart_contour_found = True

        if heart_contour_found:
            # Create mask and extract pixel values for this slice
            rr, cc = polygon(np.array([p[0] for p in img.contour_points]), np.array([p[1] for p in  img.contour_points]), img.pixel_array.shape)
            mask = np.zeros_like(img.pixel_array, dtype=bool)
            mask[rr, cc] = True
            roi_pixels = img.pixel_array[mask]

            if roi_pixels.size > 0:
                # Get the rescale slope and intercept from the CT image metadata
                rescale_slope = img.RescaleSlope
                rescale_intercept = img.RescaleIntercept

                # Convert pixel values to Hounsfield Units using the formula: HU = pixel_value * rescale_slope + rescale_intercept
                hounsfield_units = roi_pixels * rescale_slope + rescale_intercept

                # Calculate the mean Hounsfield Unit value for the heart on this CT slice
                mean_hu = np.mean(hounsfield_units)

                # Display the CT image along with the heart contour and the mean Hounsfield Unit value
                fig, ax = plt.subplots()
                ax.imshow(img.pixel_array, cmap='gray')
                ax.set_title(f"Mean HU: {mean_hu}")
                ax.plot([p[1] for p in img.contour_points], [p[0] for p in img.contour_points], linewidth=2, color='r')
                plt.show()
            else:
                print(f"No valid pixel values found for slice {img.SOPInstanceUID}")
        else:
            print(f"No heart contour found for slice {img.SOPInstanceUID}")

    except Exception as e:
        print(f"Error processing slice {img.SOPInstanceUID}: {str(e)}")

解决方案 #

要完成HU值修改，需在现有代码基础上添加像素值转换和文件保存逻辑，以下是修改后的完整代码：

import os
import pydicom
import numpy as np
from skimage.draw import polygon
import matplotlib.pyplot as plt

# 配置路径
ct_directory = r"D:\ct directory"
rt_struct_path = r"D:\RTS2.dcm"
# 新增：修改后DICOM的保存目录，需提前创建
output_directory = r"D:\modified_ct"
os.makedirs(output_directory, exist_ok=True)

# 加载CT影像
ct_files = [f for f in os.listdir(ct_directory) if f.endswith('.dcm')]
ct_images = [pydicom.dcmread(os.path.join(ct_directory, f)) for f in ct_files]
ct_images.sort(key=lambda x: float(x.ImagePositionPatient[2]))

# 加载RT结构文件
rt_struct = pydicom.dcmread(rt_struct_path)

# 定位心脏ROI编号
heart_roi_number = None
for roi in rt_struct.StructureSetROISequence:
    if 'heart' in roi.ROIName.lower():
        heart_roi_number = roi.ROINumber
        break

if heart_roi_number is None:
    raise ValueError("未找到心脏ROI")

# 处理每个CT切片
for idx, img in enumerate(ct_images):
    try:
        contour_points = []
        heart_contour_found = False

        # 匹配当前切片对应的心脏轮廓
        for roi_contour in rt_struct.ROIContourSequence:
            if roi_contour.ReferencedROINumber == heart_roi_number:
                for contour in roi_contour.ContourSequence:
                    contour_data = np.array(contour.ContourData).reshape((-1, 3))
                    z_value = contour_data[0, 2]  # 取轮廓的Z坐标（同一切片的Z值一致）
                    if abs(img.ImagePositionPatient[2] - z_value) < 1e-2:
                        # 转换轮廓点到像素坐标
                        x_coords = contour_data[:, 0]
                        y_coords = contour_data[:, 1]
                        cols = np.round((x_coords - img.ImagePositionPatient[0]) / img.PixelSpacing[0]).astype(int)
                        rows = np.round((y_coords - img.ImagePositionPatient[1]) / img.PixelSpacing[1]).astype(int)
                        contour_points.extend(list(zip(rows, cols)))
                        heart_contour_found = True

        if heart_contour_found:
            # 创建心脏区域掩码
            rr, cc = polygon(np.array([p[0] for p in contour_points]), 
                             np.array([p[1] for p in contour_points]), 
                             img.pixel_array.shape)
            mask = np.zeros_like(img.pixel_array, dtype=bool)
            mask[rr, cc] = True

            if np.any(mask):
                # 计算目标HU对应的原始像素值
                rescale_slope = img.RescaleSlope
                rescale_intercept = img.RescaleIntercept
                target_pixel_value = (100 - rescale_intercept) / rescale_slope

                # 转换为当前像素数组的 dtype（避免类型不匹配）
                target_pixel_value = np.array(target_pixel_value, dtype=img.pixel_array.dtype)

                # 修改心脏区域的像素值
                img.pixel_array[mask] = target_pixel_value

                # 保存修改后的DICOM文件
                output_path = os.path.join(output_directory, f"modified_{ct_files[idx]}")
                pydicom.dcmwrite(output_path, img, write_like_original=True)
                print(f"已修改并保存切片: {output_path}")

                # 可选：展示修改后的效果
                fig, ax = plt.subplots()
                ax.imshow(img.pixel_array, cmap='gray')
                ax.set_title("修改后心脏区域HU=100")
                ax.plot([p[1] for p in contour_points], [p[0] for p in contour_points], linewidth=2, color='r')
                plt.show()
            else:
                print(f"切片 {img.SOPInstanceUID} 无有效心脏像素")
        else:
            print(f"切片 {img.SOPInstanceUID} 未找到心脏轮廓")

    except Exception as e:
        print(f"处理切片 {img.SOPInstanceUID} 出错: {str(e)}")

关键修改说明 #

1. 新增输出目录配置：指定修改后DICOM文件的保存路径，避免覆盖原始数据
2. 优化轮廓匹配逻辑：直接取轮廓的Z坐标判断切片匹配，无需遍历每个轮廓点
3. 计算目标像素值：通过HU公式逆向推导，得到对应100HU的原始像素值，并转换为匹配的数组类型
4. 替换ROI像素：利用掩码直接修改pixel_array中对应的区域
5. 保存修改后的文件：使用pydicom.dcmwrite保存，write_like_original=True确保元数据与原始文件一致

内容的提问来源于stack exchange，提问作者aseman