问题背景：在医疗手术机器人穿刺手术中，识别显影点在空间中的坐标是一个基础问题。本篇博客简单记录在处理显影点识别时的一些思路。
目前看过两种显影球的位置：
1、显影球距离身体有一定距离
2、显影球只会出现在地板上
这里详细介绍如何处理第二种所示的情形。
显影球出现在底板时，和身体相距很近，甚至有“嵌入”的情况，如下图：
在处理第一种情况时，可以用一个滑窗（比如，滑窗大小50x50，小球尺寸固定，大概占11x11，不同扫描效果有差别）去直接扫描整幅图，因为小球距离身体足够远，我们完全可以根据一些条件判断小球是否存在与滑窗。例如，在这个滑窗内，归一化后如果像素值大于45的像素个数多余45个且小于100个，就初步认为这个滑窗内有小球存在。注：这些参数是根据数据集调整的。
但是，在第二种情况下，就不能用原来的大滑窗（比如还是50x50）去扫描。
原本，我们要定位小球在空间中的位置，想法是找出小球在轴冠矢三面哪些影像出现。将这些小球出现的轴冠矢三面的影像合起来便会形成对应的“小的”立体空间，在这些立体空间中再定位小球的空间位置。
所以，我们要求出小球在轴冠矢三面哪些影像出现。
一般而言是分别对轴冠矢三面的影像进行扫描，依次判断。
但是根据医学影像建模的特点，在LPS坐标系中，如果我们能判断小球在图片中的坐标范围，就能直接推出此小球在冠面、矢面的影像顺序。例如，某一小球在轴面，纵坐标范围[h1,h2]，横坐标范围[w1,w2]。轴面影像原本尺寸为(Z,H,W).
此小球在冠面出现的影像初始位置为：
同理，在冠面出现的最后一个影像位置为：
故，此小球在冠面出现的影像范围为[h1,h2].
同理，此小球在矢面出现的影像范围为[w1,w2]。
所以，我们只需要”准确“求出小球在轴面的坐标范围。
尝试了几种预处理方法：
ex1：对原始的轴面影像进行”隔层相减“，这种方法在第一种情况预处理时效果比较好，因为一个小球在影像中是连续出现的，比如出现在第i张~第i+6张，同时根据扫描的特性，中间位置的影像上小球是最亮的，两边缘位置的影像上小球是最暗的，那么如果我们用第i+3张减去第i-1张。在结果影像上，小球会很亮。那身体部分会怎么样呢？因为身体部分在所有影像中都是连续的，不会像小球这样在它的位置上出现“突变”。所以“隔层相减”之后，身体部分的大量像素点急剧衰减，这样的“小球部分增强，身体部分减弱”是我们想要的结果。
但这个中方法不适合第二种情况，测试时发现如果此方法应用于第二种情况，结果对“间隔数”十分敏感。还是以上面的例子说明，第i+3张小球最亮，我们减去第i-1张，在第一种情况中，第i-1张小球位置是纯黑的；但在第一种情况下，对于“部分嵌入身体”的小球，第i-1张，原本的小球区域很可能是身体部分。
ex2：对原始图像进行增强，让身体部分和小球之间的界限更加明显，然后使用霍夫圆进行检测。经过处理，将界限变得明显是可以做到的，但是霍夫圆检测不适合在这种情形中检测圆形，我一开始调参很多次，但总不能找到一组合适的参数适配手头的几个数据集。最后放弃这种思路，领导说霍夫圆检测时要求小球直径不能太小，一般需要先将原图放大，且对像素很敏感。
final：还是采取滑窗
因为要在轴面比较准确的定位小球位置。又了解到小球至多只会嵌入身体一半，不会完全嵌入。所以采用“倒T形”滑窗：
滑动时只用主窗口在图像中滑动。
具体细节如下：（代码统一放到文章末尾）
1、重采样是必须的，因为不同的数据集，虽然小球的尺寸是固定的，但轴面的图像大小不同，比如，小球在720x720大小的图像中像素数量为11x11；但在300*300大小的图像中像素数量大概为5*5。如果这样的话，小球就没有一个统一的大小标准。所以我们先对轴面图像进行重采样，我这里是把每张图像采样为700x700，小球的大小约为11x11.
2、因为小球只会出现在底部，所以我们只需要对图像底部进行扫描。这里，根据手中数据集的特点，我取图像img范围是img[520:670]。这里坐标是(y,x).
3、之后就是扫描函数：（这是核心函数）
- 因为轴面开始的某些影响还没有身体部位，但是有某些不明显的伪影，归一化之后有可能被误判为小球，所以如果原图最大像素值小于100，直接pass。如下图：
- 接下来对原图进行归一化到（0，255）
- 用7*7大小的卷积核对原图进行开运算，对于受伤数据集，经过此操作小球会变成“方形”，如下图：
- 用11*11大小的窗口去扫描图像
- 因为小球出现在底部，所以我们的扫描方式：按列进行，每一列从下向上。
- 对于每一列，遇到第一个“满的”窗口为止。“满的”窗口就是窗口内像素值大于29的个数多余105个。因为第一个“满的”窗口要么是身体部分，要么是小球，每一列最多只有一个小球，无论这个窗口是不是窗口，我们都不必要扫下去了。如果有“满的”窗口，我们就对他进行判断。否则，扫到这一列最上方才停止。
- 判断逻辑，在满足“满的”窗口前题下：
- 如果窗口最大像素值小于80，pass
- 当前窗口左上角坐标为(top,left),求出当前窗口的像素均值mean，定义thresh=mean-30
- 在主窗口左右各取一个小窗口，左边的窗口左上角坐标为：(top+5,left-11)，大小为6*11；右边的窗口左上角坐标为（top+5,left+11)，大小为6*11.
- 对于左右这两个辅助窗口，如果其中像素值大于thresh的个数超过3个，pass
- 左右两边窗口如此设置是考虑到小球最多嵌入一半
- 以上条件都满足，则该主窗口内是小球。记录下满足条件的主窗口的左上角坐标top,left，那么这个小球的纵坐标范围就是(top,top+11);横坐标范围是(left,left+11).据此能推算出对应的冠状面、矢状面的位置。（实际情况会扩大一下范围）
- 行列步长都设置的为1.这样很慢，加速方式采用的积分图，大意便是利用掩模图，对于点（i，j），计算出他作为左上角时大小为11*11的窗口内大于29的值个数，最大像素值，记录。可以结合代码看。
通过以上思路，便可以计算出小球在轴冠矢三面出现的影像位置。对于100*700*700的数据集，时间不超过2s。
具体代码如下：
code 1，重采样函数：
def resample_dicom_series(dicom_dir, target_size_xy=(700, 700)):
# 1. 读取 DICOM 序列
reader = sitk.ImageSeriesReader()
dicom_names = reader.GetGDCMSeriesFileNames(dicom_dir)
reader.SetFileNames(dicom_names)
image = reader.Execute() # 3D 图像
original_size = image.GetSize() # (x, y, z)
original_spacing = image.GetSpacing() # (x, y, z)
# print(f"原始尺寸: {original_size}")
# print(f"原始spacing: {original_spacing}")
# 2. 设定新尺寸
new_size = [
target_size_xy[0],
target_size_xy[1],
original_size[2] # Z轴不变
]
# 3. 计算新的 spacing，保证物理尺寸不变
new_spacing = [
original_spacing[0] * (original_size[0] / new_size[0]),
original_spacing[1] * (original_size[1] / new_size[1]),
original_spacing[2] # Z轴 spacing 不变
]
# print(f"新尺寸: {new_size}")
# print(f"新 spacing: {new_spacing}")
# 4. 重采样
resampler = sitk.ResampleImageFilter()
resampler.SetSize(new_size)
resampler.SetOutputSpacing(new_spacing)
resampler.SetOutputOrigin(image.GetOrigin())
resampler.SetOutputDirection(image.GetDirection())
resampler.SetInterpolator(sitk.sitkLinear)
resampled_image = resampler.Execute(image)
resampled_image = sitk.Cast(resampled_image, sitk.sitkUInt16) # 转换为uint16
# 可选：保存为 NumPy 数组
np_array = sitk.GetArrayFromImage(resampled_image) # shape: [z, y, x]
size = image.GetSize()
return (size[1], size[0]), np_array
扫描函数：（核心代码）
def sliding_sum(mask, k):
"""mask 为 uint8/uint16/float32 图像，返回每个 k×k 窗口的像素和，左上角对齐"""
intg = cv2.integral(mask, sdepth=cv2.CV_32S) # (h+1,w+1)
# h, w = mask.shape
out = (
intg[k:, k:] - intg[:-k, k:] -
intg[k:, :-k] + intg[:-k, :-k]
)
return out # 形状为 (h-k+1, w-k+1)
def sliding_max(img, k):
kernel = np.ones((k, k), np.uint8)
return cv2.dilate(img, kernel, anchor=(0, 0))[0:img.shape[0] - k + 1, 0:img.shape[1]
- k + 1]
def scan_image_vec(vol, k=11):
h, w = vol[0].shape
H, W = h - k + 1, w - k + 1 # 有效窗口行/列数
coords = [] # 存储结果坐标
for i in range(0, len(vol)):
img = vol[i]
if np.max(img) < 100: # 如果原图几乎是黑色的，直接排除
continue
# 归一化
norm_img = cv2.normalize(img, None, 0, 255, cv2.NORM_MINMAX)
# 开运算
kernel = cv2.getStructuringElement(cv2.MORPH_RECT, (7, 7))
img = cv2.morphologyEx(norm_img, cv2.MORPH_OPEN, kernel)
cnt19 = sliding_sum((img > 29).astype(np.uint8), k)
win_max = sliding_max(img, k)
win_sum = sliding_sum(img.astype(np.uint8), k)
win_mean = win_sum / (k * k) - 30
full = (cnt19 >= 105) # 记录满的窗口
bigger_80 = win_max >= 80 # 记录最大值大于80的窗口
full_flip = full[::-1, :] # 上下对调
first_true_idx_from_bottom = np.argmax(full_flip, axis=0)
col = k - 1
while col < W:
# 自下而上找本列第一处候选
# rows = np.flatnonzero(full[::-1, col]) # 求出full矩阵每一列从下向上第一个“满的”矩阵位置
row = first_true_idx_from_bottom[col]
if row == 0: # 如果此列没有满的矩阵，直接跳过
col += 1
continue
top = H - 1 - row # 转回正向行号
if not bigger_80[top, col]: # 如果这个滑窗正好最大值<80,pass
col += 1
continue
thresh = win_mean[top, col]
next_top = top + 5
# ---------- 左右 6×11 判定 ----------
ok = False
# 左侧 右侧
l0 = col - k
r0 = col + k
if l0 >= 0 and r0 < w:
left_win = img[next_top:next_top + 6, l0:l0 + k]
right_win = img[next_top:next_top + 6, r0:r0 + k]
if (left_win > thresh).sum() <= 3 and (right_win > thresh).sum() <= 3:
ok = True
# -------------------------------------
if ok: # 真·小球
coords.append((i, top, col))
col += 20 # 两个相邻小球列间隔不会小于20
else:
col += 1
return coords
下面的函数是推算冠矢面位置的代码：
def select_idx(valid_coords):
cur = [(p[1],p[2]) for p in valid_coords]
cur.sort(key=lambda it: (it[1], it[0]))
ans = [cur[0]]
for i in range(1, len(cur)):
if abs(cur[i][0] - cur[i-1][0]) <= 2 and abs(cur[i][1] - cur[i-1][1]) <= 2:
continue
else:
ans.append(cur[i])
return ans
扫一扫
37
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
