python实现血管增强算法

血管增强算法是医学图像处理中的重要技术，可以提高血管的对比度，辅助医生进行诊断。以下是一个基于 Python 的血管增强算法实现：

1.导入相关库

import cv2
import numpy as np
from matplotlib import pyplot as plt

2.读取原始图像

img = cv2.imread('original_image.png',0)

3.预处理

对原始图像进行预处理，去除噪声和平滑化。

img = cv2.medianBlur(img,5)
img = cv2.GaussianBlur(img,(5,5),0)

4.计算梯度幅值和方向

计算图像的梯度幅值和方向，用于后续处理。

sobelx = cv2.Sobel(img,cv2.CV_64F,1,0,ksize=3)
sobely = cv2.Sobel(img,cv2.CV_64F,0,1,ksize=3)
mag = np.sqrt(sobelx**2 + sobely**2)
theta = np.arctan2(sobely, sobelx)

5.非极大值抑制

在梯度方向上进行非极大值抑制，保留局部最大值点。

rows,cols = img.shape
for i in range(1,rows-1):
    for j in range(1,cols-1):
        if (theta[i,j] >= -np.pi/8 and theta[i,j] < np.pi/8) or (theta[i,j] >= 7*np.pi/8 and theta[i,j] <= np.pi):
            if (mag[i,j] >= mag[i,j-1]) and (mag[i,j] >= mag[i,j+1]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0
        elif (theta[i,j] >= np.pi/8 and theta[i,j] < 3*np.pi/8):
            if (mag[i,j] >= mag[i-1,j+1]) and (mag[i,j] >= mag[i+1,j-1]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0
        elif (theta[i,j] >= 3*np.pi/8 and theta[i,j] < 5*np.pi/8):
            if (mag[i,j] >= mag[i-1,j]) and (mag[i,j] >= mag[i+1,j]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0
        else:
            if (mag[i,j] >= mag[i-1,j-1]) and (mag[i,j] >= mag[i+1,j+1]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0

6.双阈值处理

根据梯度幅值大小进行双阈值处理，将图像分成强边缘、弱边缘和非边缘三类。

high_threshold, thresh_im = cv2.threshold(mag, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
low_threshold = 0.5*high_threshold
rows,cols = img.shape
for i in range(rows):
    for j in range(cols):
        if mag[i,j] >= high_threshold:
            mag[i,j] = 255
        elif mag[i,j] >= low_threshold and mag[i,j] < high_threshold:
            mag[i,j] = 100
        else:
            mag[i,j] = 0

7.连接边缘

对强边缘进行连接，得到完整的血管。

for i in range(1,rows-1):
    for j in range(1,cols-1):
        if mag[i,j] == 255:
            if mag[i-1,j-1] == 100:
                mag[i-1,j-1] = 255
            if mag[i-1,j] == 100:
                mag[i-1,j] = 255
            if mag[i-1,j+1] == 100:
                mag[i-1,j+1] = 255
            if mag[i,j-1] == 100:
                mag[i,j-1] = 255
            if mag[i,j+1] == 100:
                mag[i,j+1] = 255
            if mag[i+1,j-1] == 100:
                mag[i+1,j-1] = 255
            if mag[i+1,j] == 100:
                mag[i+1,j] = 255
            if mag[i+1,j+1] == 100:
                mag[i+1,j+1] = 255

8.显示结果

plt.subplot(121),plt.imshow(img,cmap = 'gray')
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(mag,cmap = 'gray')
plt.title('Enhanced Image'), plt.xticks([]), plt.yticks([])
plt.show()

完整代码如下：

import cv2
import numpy as np
from matplotlib import pyplot as plt

img = cv2.imread('original_image.png',0)

img = cv2.medianBlur(img,5)
img = cv2.GaussianBlur(img,(5,5),0)

sobelx = cv2.Sobel(img,cv2.CV_64F,1,0,ksize=3)
sobely = cv2.Sobel(img,cv2.CV_64F,0,1,ksize=3)
mag = np.sqrt(sobelx**2 + sobely**2)
theta = np.arctan2(sobely, sobelx)

rows,cols = img.shape
for i in range(1,rows-1):
    for j in range(1,cols-1):
        if (theta[i,j] >= -np.pi/8 and theta[i,j] < np.pi/8) or (theta[i,j] >= 7*np.pi/8 and theta[i,j] <= np.pi):
            if (mag[i,j] >= mag[i,j-1]) and (mag[i,j] >= mag[i,j+1]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0
        elif (theta[i,j] >= np.pi/8 and theta[i,j] < 3*np.pi/8):
            if (mag[i,j] >= mag[i-1,j+1]) and (mag[i,j] >= mag[i+1,j-1]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0
        elif (theta[i,j] >= 3*np.pi/8 and theta[i,j] < 5*np.pi/8):
            if (mag[i,j] >= mag[i-1,j]) and (mag[i,j] >= mag[i+1,j]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0
        else:
            if (mag[i,j] >= mag[i-1,j-1]) and (mag[i,j] >= mag[i+1,j+1]):
                mag[i,j] = mag[i,j]
            else:
                mag[i,j] = 0

high_threshold, thresh_im = cv2.threshold(mag, 0, 255, cv2.THRESH_BINARY + cv2.THRESH_OTSU)
low_threshold = 0.5*high_threshold
rows,cols = img.shape
for i in range(rows):
    for j in range(cols):
        if mag[i,j] >= high_threshold:
            mag[i,j] = 255
        elif mag[i,j] >= low_threshold and mag[i,j] < high_threshold:
            mag[i,j] = 100
        else:
            mag[i,j] = 0

for i in range(1,rows-1):
    for j in range(1,cols-1):
        if mag[i,j] == 255:
            if mag[i-1,j-1] == 100:
                mag[i-1,j-1] = 255
            if mag[i-1,j] == 100:
                mag[i-1,j] = 255
            if mag[i-1,j+1] == 100:
                mag[i-1,j+1] = 255
            if mag[i,j-1] == 100:
                mag[i,j-1] = 255
            if mag[i,j+1] == 100:
                mag[i,j+1] = 255
            if mag[i+1,j-1] == 100:
                mag[i+1,j-1] = 255
            if mag[i+1,j] == 100:
                mag[i+1,j] = 255
            if mag[i+1,j+1] == 100:
                mag[i+1,j+1] = 255

plt.subplot(121),plt.imshow(img,cmap = 'gray')
plt.title('Original Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(mag,cmap = 'gray')
plt.title('Enhanced Image'), plt.xticks([]), plt.yticks([])
plt.show()