如何定位DICOM图像中最大分割区域的中心？

我按照一份DICOM图像分析指南的代码完成步骤3后，得到了切片图像的分割区域（见下图）。现在需要找到右侧图像中白色大区域的中心。

我不确定有没有直接定位白色区域中心的方法，若有这类方法会很有帮助。我的初步思路是通过测量区域的宽高，在图像上画两条垂直线来确定中心，但不确定这个方法的效率，也不清楚具体实现方式。

任何能测量分割区域中心的方法我都很感激，谢谢。

以下是我生成分割图像的代码：

import keras, os, copy, pydicom, cv2, glob
import matplotlib.pyplot as plt
import numpy as np
import nibabel as nib
import scipy.ndimage
import pandas as pd
from skimage import filters
from math import *
from skimage import measure, morphology
from skimage.morphology import ball, binary_closing
from skimage.measure import label, regionprops
from functools import reduce
from scipy.linalg import norm

path = "filepath" #for me, the filepath is to a DICOM .dcm file

def load_scan(path):
    slices = [pydicom.dcmread(path + '/' + s) for s in               
              os.listdir(path)]
    slices = [s for s in slices if 'SliceLocation' in s]
    slices.sort(key = lambda x: int(x.InstanceNumber))
    try:
        slice_thickness = np.abs(slices[0].ImagePositionPatient[2] -   
                          slices[1].ImagePositionPatient[2])
    except:
        slice_thickness = np.abs(slices[0].SliceLocation - 
                      slices[1].SliceLocation)
    for s in slices:
        s.SliceThickness = slice_thickness
    return slices

def get_pixels_hu(scans):
    image = np.stack([s.pixel_array for s in scans])
    image = image.astype(np.int16)
    # Set outside-of-scan pixels to 0
    # The intercept is usually -1024, so air is approximately 0
    image[image == -2000] = 0
    
    # Convert to Hounsfield units (HU)
    intercept = scans[0].RescaleIntercept
    slope = scans[0].RescaleSlope
    
    if slope != 1:
        image = slope * image.astype(np.float64)
        image = image.astype(np.int16)
        
    image += np.int16(intercept)
    
    return np.array(image, dtype=np.int16)

patient_dicom = load_scan(path)
patient_pixels = get_pixels_hu(patient_dicom)

def largest_label_volume(im, bg=-1):
    vals, counts = np.unique(im, return_counts=True)
    counts = counts[vals != bg]
    vals = vals[vals != bg]
    if len(counts) > 0:
        return vals[np.argmax(counts)]
    else:
        return None
    
def segment_lung_mask(image, fill_lung_structures=True):
    # not actually binary, but 1 and 2. 
    # 0 is treated as background, which we do not want
    binary_image = np.array(image >= -800, dtype = np.int8) + 1
    labels = measure.label(binary_image)
 
    # Pick the pixel in the very corner to determine which label is air.
    # Improvement: Pick multiple background labels from around the patient
    # More resistant to “trays” on which the patient lays cutting the air around the person in half 
    background_label = labels[0,0,0]
 
    # Fill the air around the person
    binary_image[background_label == labels] = 1
 
    # Method of filling the lung structures (that is superior to 
    # something like morphological closing)
    if fill_lung_structures:
        # For every slice we determine the largest solid structure
        for i, axial_slice in enumerate(binary_image):
            axial_slice = axial_slice - 1
            labeling = measure.label(axial_slice)
            l_max = largest_label_volume(labeling, bg=0)
 
            if l_max is not None: #This slice contains some lung
                binary_image[i][labeling != l_max] = 1
            
    binary_image -= 1 #Make the image actual binary
    binary_image = 1-binary_image # Invert it, lungs are now 1
 
    # Remove other air pockets inside body
    labels = measure.label(binary_image, background=0)
    l_max = largest_label_volume(labels, bg=0)
    if l_max is not None: # There are air pockets
        binary_image[labels != l_max] = 0
 
    return binary_image

# get masks 
segmented_lungs = segment_lung_mask(patient_pixels,    
                    fill_lung_structures = False)
segmented_lungs_fill = segment_lung_mask(patient_pixels,     
                    fill_lung_structures = True)
#internal_structures = segmented_lungs_fill - segmented_lungs

# isolate lung from chest
copied_pixels = copy.deepcopy(patient_pixels)
for i, mask in enumerate(segmented_lungs_fill): 
    get_high_vals = mask == 0
    copied_pixels[i][get_high_vals] = 0
seg_lung_pixels = copied_pixels

# sanity check
f, ax = plt.subplots(1,2, figsize=(10,6))
ax[0].imshow(patient_pixels[60], cmap=plt.cm.bone)
ax[0].axis(False)
ax[0].set_title('Original')
ax[1].imshow(seg_lung_pixels[60], cmap=plt.cm.bone)
ax[1].axis(False)
ax[1].set_title('Segmented')
plt.show()

要定位分割后白色区域的中心，有两种高效且易实现的方法，适配医学图像的不规则区域特性：

方法一：计算区域质心（精准贴合形状） #

质心是区域的加权中心，能准确反映不规则区域的核心位置，适合肺部这类非规则分割区域。结合现有代码，实现步骤如下：

1. 提取目标切片的分割掩码：

# 取第60层的二值分割掩码（1代表肺部区域）
slice_mask = segmented_lungs_fill[60]

2. 分析区域属性并获取质心：

# 标记区域并提取属性
labeled_mask = measure.label(slice_mask)
region_props = regionprops(labeled_mask)

# 筛选面积最大的区域（对应图像中的白色大区域）
largest_region = max(region_props, key=lambda prop: prop.area)
# 注意：skimage返回的质心格式为(y, x)，对应图像的行、列索引
centroid_y, centroid_x = largest_region.centroid

3. 在图像上可视化中心：

slice_image = seg_lung_pixels[60]

plt.figure(figsize=(8,8))
plt.imshow(slice_image, cmap=plt.cm.bone)
# 绘制质心标记
plt.scatter(centroid_x, centroid_y, color='red', s=50, marker='x', label='质心')
# 添加辅助线
plt.axvline(x=centroid_x, color='blue', linestyle='--')
plt.axhline(y=centroid_y, color='blue', linestyle='--')
plt.legend()
plt.axis(False)
plt.title('分割区域质心')
plt.show()

方法二：计算几何中心（快速定位） #

如果只需要大致的中心位置，可以基于区域的最小外接矩形计算几何中心，速度更快：

# 从最大区域获取边界框坐标
min_row, min_col, max_row, max_col = largest_region.bbox
# 计算几何中心
center_y = (min_row + max_row) / 2
center_x = (min_col + max_col) / 2

# 可视化几何中心
plt.figure(figsize=(8,8))
plt.imshow(slice_image, cmap=plt.cm.bone)
plt.scatter(center_x, center_y, color='green', s=50, marker='o', label='几何中心')
plt.axvline(x=center_x, color='purple', linestyle='--')
plt.axhline(y=center_y, color='purple', linestyle='--')
plt.legend()
plt.axis(False)
plt.title('分割区域几何中心')
plt.show()

注意事项 #

• 质心考虑了区域内像素的分布，精度更高，适合医学图像分析场景；
• 几何中心计算逻辑简单，适合快速获取大致位置；
• 两种方法均可直接集成到现有代码的sanity check代码之后，无需修改核心分割逻辑。

内容的提问来源于stack exchange，提问作者Just A Mathematician