一种模拟金属伪影的方法代码

贴张演示图

 1、读取图片，我这里是dcm

dcm_path = ""
dcm = pydicom.read_file(dcm_path)
gt = dcm.pixel_array

2、对每一个slice进行归一化，加金属伪影，反归一化

withMA = np.zeros_like(gt)
for i in range(gt.shape[0]):
    print(i)
    ori = normalized(gt[i,:,:])
    with_mask=Ma(ori) * (4000) - 1000
    withMA[i] =with_mask

3、简单的归一化

def normalized(X):
    maxX = 3000
    minX = -1000
    X=np.clip(X, minX, maxX)
    X = (X - minX) / (maxX - minX)
    return X

4、金属伪影，需要调整的：

vec =(-shape/2,shape/2),

x1,y1,为伪影位置，r1为伪影半径

a,c为伪影强度

def Ma(gt):
    # 建立金属掩码
    metalMask=np.zeros_like(gt)
    vec=np.arange(-135, 135)
    xx, yy=np.meshgrid(vec, vec)
    x1=random.randint(-20, 20);
    y1=-82 + np.abs(x1);
    r1=4

    metalMask[(xx - x1) ** 2 + (yy - y1) ** 2 <= r1 ** 2]=1

    # 金属和非金属部分的Radon变换
    radon_theta=np.arange(0, 180)
    R_metal=radon(metalMask, theta=radon_theta, circle=True)
    R_withMetal=radon(gt, theta=radon_theta, circle=True) + R_metal
    # test_image=iradon(R_withMetal, theta=radon_theta, circle=True)

    # 金属投影区域进行模拟变换
    prj_metalSim=R_withMetal.copy()
    a=6
    c=0.02

    # 假设 x 和 h 是根据某种方式计算得到的
    x=np.linspace(a - 5, a + 5, 10)  # 举例的 x 值范围
    h=x - (2 * c * (x - a) + 1 - np.sqrt(np.maximum(4 * c * (x - a) + 1, 0))) / (2 * c) * np.sign(x - a)

    # 参数 c 控制变换曲线的宽度，即金属伪影的过渡区域的宽度。
    # 参数 a 是变换曲线中心的位置。通过调整这些参数的值，可以改变 h 函数的形状，进而影响金属伪影的外观。
    # 创建插值函数，使用 fill_value='extrapolate' 来允许超出范围的插值
    h_func=interp1d(x, h, kind='cubic', fill_value="extrapolate")

    for i in range(prj_metalSim.shape[0]):
        prj_metalSim[i, :]=h_func(prj_metalSim[i, :])

    metal_image=iradon(prj_metalSim, theta=radon_theta, circle=True)
    return metal_image

5、展示

plt.figure();plt.imshow(withMA[0], cmap='gray');plt.title('模拟金属伪影的CT图')
plt.show()