2

我目前正在尝试在 Python 中实现这种用于体绘制的算法，并且在概念上对他们生成 LH 直方图的方法感到困惑（参见第 3.1 节，第 4 页）。

我有一个 DICOM 图像的 3D 堆栈，并用它计算了它的梯度幅度和 2 个相应的方位角和仰角（我在这里找到了），以及找到二阶导数。

现在，该算法要求我遍历一组体素，并“通过在两个方向上积分梯度场来跟踪路径……使用具有一个体素积分步骤的二阶龙格-库塔方法”。

我不明白的是如何使用我计算出的 2 个角度来整合所述方向的梯度场。我知道您可以使用三线性插值来获得中间体素值，但我不明白如何使用我拥有的角度获得我想要的体素坐标。

换句话说，我从给定的体素位置开始，并希望在为该体素计算的 2 个角度的方向上采取 1 个体素步骤（一个在 xy 方向，另一个在 z 方向）。我将如何在这两个角度采取这一步并检索新的 (x, y, z) 体素坐标？

提前道歉，因为我在 Calc II/III 方面有非常基本的背景，所以矢量场/3D 空间的可视化对我来说还是有点粗糙。

创建 DICOM 图像的 3D 堆栈：

def collect_data(data_path):
    print "collecting data"
    files = []  # create an empty list
    for dirName, subdirList, fileList in os.walk(data_path):
        for filename in fileList:
            if ".dcm" in filename:
                files.append(os.path.join(dirName,filename))
    # Get reference file
    ref = dicom.read_file(files[0])

    # Load dimensions based on the number of rows, columns, and slices (along the Z axis)
    pixel_dims = (int(ref.Rows), int(ref.Columns), len(files))

    # Load spacing values (in mm)
    pixel_spacings = (float(ref.PixelSpacing[0]), float(ref.PixelSpacing[1]), float(ref.SliceThickness))

    x = np.arange(0.0, (pixel_dims[0]+1)*pixel_spacings[0], pixel_spacings[0])
    y = np.arange(0.0, (pixel_dims[1]+1)*pixel_spacings[1], pixel_spacings[1])
    z = np.arange(0.0, (pixel_dims[2]+1)*pixel_spacings[2], pixel_spacings[2])

    # Row and column directional cosines
    orientation = ref.ImageOrientationPatient

    # This will become the intensity values
    dcm = np.zeros(pixel_dims, dtype=ref.pixel_array.dtype)

    origins = []

    # loop through all the DICOM files
    for filename in files:
        # read the file
        ds = dicom.read_file(filename)
        #get pixel spacing and origin information
        origins.append(ds.ImagePositionPatient) #[0,0,0] coordinates in real 3D space (in mm)

        # store the raw image data
        dcm[:, :, files.index(filename)] = ds.pixel_array
    return dcm, origins, pixel_spacings, orientation

计算梯度幅度：

def calculate_gradient_magnitude(dcm):
    print "calculating gradient magnitude"
    gradient_magnitude = []
    gradient_direction = []

    gradx = np.zeros(dcm.shape)
    sobel(dcm,0,gradx)
    grady = np.zeros(dcm.shape)
    sobel(dcm,1,grady)
    gradz = np.zeros(dcm.shape)
    sobel(dcm,2,gradz)

    gradient = np.sqrt(gradx**2 + grady**2 + gradz**2)

    azimuthal = np.arctan2(grady, gradx)
    elevation = np.arctan(gradz,gradient)

    azimuthal = np.degrees(azimuthal)
    elevation = np.degrees(elevation)

    return gradient, azimuthal, elevation

转换为患者坐标系以获得实际体素位置：

def get_patient_position(dcm, origins, pixel_spacing, orientation):
    """
        Image Space --> Anatomical (Patient) Space is an affine transformation
        using the Image Orientation (Patient), Image Position (Patient), and
        Pixel Spacing properties from the DICOM header
    """
    print "getting patient coordinates"

    world_coordinates = np.empty((dcm.shape[0], dcm.shape[1],dcm.shape[2], 3))
    affine_matrix = np.zeros((4,4), dtype=np.float32)

    rows = dcm.shape[0]
    cols = dcm.shape[1]
    num_slices = dcm.shape[2]

    image_orientation_x = np.array([ orientation[0], orientation[1], orientation[2]  ]).reshape(3,1)
    image_orientation_y = np.array([ orientation[3], orientation[4], orientation[5]  ]).reshape(3,1)
    pixel_spacing_x = pixel_spacing[0]

    # Construct affine matrix
    # Method from:
    # http://nipy.org/nibabel/dicom/dicom_orientation.html
    T_1 = origins[0]
    T_n = origins[num_slices-1]


    affine_matrix[0,0] = image_orientation_y[0] * pixel_spacing[0]
    affine_matrix[0,1] = image_orientation_x[0] * pixel_spacing[1]
    affine_matrix[0,3] = T_1[0]
    affine_matrix[1,0] = image_orientation_y[1] * pixel_spacing[0]
    affine_matrix[1,1] = image_orientation_x[1] * pixel_spacing[1]
    affine_matrix[1,3] = T_1[1]
    affine_matrix[2,0] = image_orientation_y[2] * pixel_spacing[0]
    affine_matrix[2,1] = image_orientation_x[2] * pixel_spacing[1]
    affine_matrix[2,3] = T_1[2]
    affine_matrix[3,3] = 1

    k1 = (T_1[0] - T_n[0])/ (1 - num_slices)
    k2 = (T_1[1] - T_n[1])/ (1 - num_slices)
    k3 = (T_1[2] - T_n[2])/ (1 - num_slices)

    affine_matrix[:3, 2] = np.array([k1,k2,k3])

    for z in range(num_slices):
        for r in range(rows):
            for c in range(cols):
                vector = np.array([r, c, 0, 1]).reshape((4,1))
                result = np.matmul(affine_matrix, vector)
                result = np.delete(result, 3, axis=0)
                result = np.transpose(result)
                world_coordinates[r,c,z] = result
        # print "Finished slice ", str(z)
    # np.save('./data/saved/world_coordinates_3d.npy', str(world_coordinates))
    return world_coordinates

现在我要编写这个函数了：

def create_lh_histogram(patient_positions, dcm, magnitude, azimuthal, elevation):
    print "constructing LH histogram"
    # Get 2nd derivative
    second_derivative = gaussian_filter(magnitude, sigma=1, order=1)

    # Determine if voxels lie on boundary or not (thresholding)

    # Still have to code out: let's say the thresholded voxels are in 
    # a numpy array called voxels

    #Iterate through all thresholded voxels and integrate gradient field in
    # both directions using 2nd-order Runge-Kutta
    vox_it = voxels.nditer(voxels, flags=['multi_index'])
    while not vox_it.finished:
        # ???
4

0 回答 0