Hounsfield 单位 (DICOM) 的像素值

Question

我有一个切片卷文件 (.stl) 并将切片转换为 DICOM 文件的代码。我得到二值化图像（0 或 255），我想将它们转换为 DICOM 文件，对于 255 的像素，Hounsfield 单位为 0，对于 0 的像素，则为 -1000。这是代码的一部分那：

# Previous image (binary)
img = ndarray_to_pil(img).convert("1")  
arr = pil_to_ndarray(img) 

# Convert to Hounsfield units
arr = arr.astype(np.float32)
arr[arr == 0] = -1000
arr[arr == 255] = 0

# Create a DICOM image
# Set the image pixel data and metadata
ds.Modality = 'CT'
ds.PixelData = arr.tobytes()
ds.SamplesPerPixel = 1
ds.PhotometricInterpretation = 'MONOCHROME2'
ds.Rows, ds.Columns = arr.shape[:2]
ds.BitsAllocated = 8
ds.BitsStored = 8
ds.HighBit = 7
ds.PixelRepresentation = 0
ds.RescaleIntercept = np.min(arr)
ds.RescaleSlope = 1.0
ds.PixelSpacing = [1.0, 1.0]
ds.ImagePositionPatient = [0.0, 0.0, 0.0]
ds.ImageOrientationPatient = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]

# Set the window level and width (adjust these values as needed)
ds.WindowCenter = np.max(arr) / 2
ds.WindowWidth = np.max(arr)

# Save the DICOM file
ds.save_as(os.path.join(slice_folder, f'slice_{i}.dcm'))

这是原始的.png二值化图像

但这是我得到的结果（我用ImageJ打开它）。他们的 Haunsfield 值有奇怪的条纹，分别是 -1000、-773 和 -792：

（完整代码）：

import numpy as np
import pyvista as pv
import os
import matplotlib
import cv2
matplotlib.use('Agg') 
import matplotlib.pyplot as plt
from PIL import Image
from skimage.morphology import skeletonize, thin
from skimage.io._plugins.pil_plugin import ndarray_to_pil, pil_to_ndarray
import tkinter as tk
import pydicom
from tkinter import filedialog
import warnings

warnings.filterwarnings("ignore", category=FutureWarning)


root = tk.Tk()
root.title('Tkinter Open File Dialog')
root.resizable(False, False)
root.geometry('300x150')


stl_path = filedialog.askopenfilename(
    title="Selecciona un archivo STL",
    filetypes=(("STL files", "*.stl"), ("All files", "*.*"))
)

slice_folder = filedialog.askdirectory(
    title="Selecciona una carpeta donde guardar el stack de imágenes",
)

root.withdraw()


stl_file = pv.read(stl_path)

xmin, xmax, ymin, ymax, zmin, zmax = stl_file.bounds

num_slices = 32

slice_height = (zmax - zmin) / num_slices


# Slicing loop
for i in range(num_slices):
    z = zmin + i * slice_height
    mesh_slice = stl_file.slice_along_axis(n=1, axis='z', bounds=(0, 0, 0, 0, z, z), generate_triangles=True)
    mesh_slice = mesh_slice[0]  # Extract the first block
    x, y, z = mesh_slice.points.transpose()
    fig = plt.figure()
    plt.style.use('dark_background')
    plt.plot(x, y, 'o', markersize=12)
    plt.axis([xmin, xmax, ymin, ymax]) # rango del plano XY del .stl en el plot
    plt.axis('off')
    plt.close('all')

    fig.canvas.draw()
    img = Image.frombytes('RGB', fig.canvas.get_width_height(), fig.canvas.tostring_rgb())
    
    img = cv2.blur(np.array(img),(5,5))
    img = Image.fromarray(img.astype(np.uint8))
    
    # Skeletization
    img_gray = img.convert('L')    
    img = img_gray.resize((512, 512), resample=Image.LANCZOS)
    ret, img_binary = cv2.threshold(np.array(img), 10, 255, 0)
    img = skeletonize(img_binary)
    img = ndarray_to_pil(img).convert("1")  
    img_array = pil_to_ndarray(img)
    img_array = img_array.reshape(512, 512, 1)
    img = img_array[:,:,0]
    ret, thresh = cv2.threshold(img, 0, 255, cv2.THRESH_BINARY)
    
    
    # Find the outer contours in the binary image (using cv2.RETR_EXTERNAL)
    contours, hierarchy = cv2.findContours(thresh, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
    
    # Create a blank image with the same dimensions as the original image
    filled_img = np.zeros(img.shape[:2], dtype=np.uint8)
    
    # Fill the outer contour with white color
    cv2.drawContours(filled_img, contours, -1, 255, cv2.FILLED)

    # Find contours with hierarchy, this time use cv2.RETR_TREE
    contours, hierarchy = cv2.findContours(thresh, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
    
    # Iterate over the contours and their hierarchies
    for j, contour in enumerate(contours):
        has_grandparent = False
        has_parent = hierarchy[0][j][3] >= 0
        if has_parent:
            # Check if contour has a grandparent
            parent_idx = hierarchy[0][j][3]
            has_grandparent = hierarchy[0][parent_idx][3] >= 0
    
        # Check if the contour has no child
        if hierarchy[0][j][2] < 0 and has_grandparent:
            # If contour has no child, fill the contour with black color
            cv2.drawContours(filled_img, [contour], -1, 0, cv2.FILLED)
    
    img = ndarray_to_pil(filled_img).convert("1")  
    arr = pil_to_ndarray(img)     
    
    ### DICOM CONVERSION ###
    
    # Haunsfield units
    arr = arr.astype(np.float16)
    arr[arr == 0] = -1000
    arr[arr == 255] = 0

    
    # Create dicom file
    ds = pydicom.dataset.FileDataset(os.path.join(slice_folder, f'slice_{i}.dcm'), {}, file_meta=pydicom.dataset.FileMetaDataset())
    
    # Metadata
    ds.PatientName = 'Nombre'
    ds.PatientID = 'ID'

    # Set the image pixel data and metadata
    ds.Modality = 'CT'
    ds.PixelData = arr.tobytes()
    ds.SamplesPerPixel = 1
    ds.PhotometricInterpretation = 'MONOCHROME2'
    ds.Rows, ds.Columns = arr.shape[:2]
    ds.BitsAllocated = 8
    ds.BitsStored = 8
    ds.HighBit = 7
    ds.PixelRepresentation = 0
    ds.RescaleIntercept = np.min(arr)
    ds.RescaleSlope = 1.0
    ds.PixelSpacing = [1.0, 1.0]
    ds.ImagePositionPatient = [0.0, 0.0, 0.0]
    ds.ImageOrientationPatient = [1.0, 0.0, 0.0, 0.0, 1.0, 0.0]

    # Set the window level and width (adjust these values as needed)
    ds.WindowCenter = np.max(arr) / 2
    ds.WindowWidth = np.max(arr)

    # Save the DICOM file
    ds.save_as(os.path.join(slice_folder, f'slice_{i}.dcm'))
    
root.destroy()
print('done')