图像血管提取python 血管图像后处理技术

    终于把这个Project弄完了。我采取的是Python的轻量级的图像处理插件PIL进行包括图像增强，图像分割等方面的处理。实验室采集的手背静脉图像低对比度低分辨率噪声非常大。整个处理过程包括预处理，ROI定位，以及后续处理，实验结果还算理想。

   原始图片：

  

（此图如有版权问题请与我联系）

ROI：

最后的细化结果:

    网上有很多人采取的是Python与Opencv结合的方式进行数字图像的处理，更多的是直接使用matlab。由于我采用的是纯python语言，因此我做了很多PIL没有做到的但是Opencv做的事情。PIL是轻量级的Python Image 处理工具，用于专业的数字图像处理还是不怎么适合的，在这个过程中我花费了很多时间，也算收获不少吧。

   基于手背血管进行身份识别的预处理研究  是老师给我们的一篇参考论文，作者是天津理工大学的***，这篇论文竟然跟清华大学的***写的论文

人体手背血管图像的特征提取及匹配惊人地相似，我也无法甄别其中的真伪，但有一点可以肯定的是一方过分参考另外一方，呵呵，在中国，特别是搞学术的这个大家都懂的，这其中的原因很复杂，博客园是个和谐的地方～～～中国的学术环境挺恶劣的说。

1.直方图均衡化：

     很经典，大家熟知的图像增强方法。 

# -*- coding: utf-8 -*-

# histogram equalization 还得写均衡化程序晕死

import operator
import Image
def equalize(h):
    lut = []
for b in range(0, len(h), 256):
# 步数
        step = reduce(operator.add, h[b:b+256]) / 255
# 创建查找表
        n = 0
for i in range(256):
            lut.append(n / step)
            n = n + h[i+b]
return lut

#测试程序
if __name__ == "__main__": 
    im = Image.open('roi.jpg')
#计算查找表
    lut = equalize(im.histogram())
# 查表
    im = im.point(lut)
    im.save("roi_equalized.jpg")

2.二值化：

# -*- coding: utf-8 -*-
#Binary operation using a given threshold
#Author:[e-mail] 16时51分31秒20110516

import Image

def fxxBinary(im, threshold): #给定阈值，返回二值图像
    # convert to grey level image
    Lim = im.convert('L')
    Lim.save('roi_Level.jpg')

# setup a converting table with constant threshold
    table = []
for i in range(256):
if i < threshold:
            table.append(0)
else:
            table.append(255)

# convert to binary image by the table
    bim = Lim.point(table, '1')
return bim

#测试代码
if __name__ =='__main__':
    im = Image.open('roi_equalized.jpg')
    im.show()
    threshold = 128
    bim = fxxBinary(im, threshold)
    bim.show()
    bim.save('roi_binary.jpg')

3  Otsu阈值分割 

     matlab中的im2bw采用的就是这种阈值分割方式。借助于上面的二值化代码。

# -*- coding: utf-8 -*-
import Image,ImageEnhance,ImageFilter,ImageDraw
import fxxBinary
'''大津法由大津于1979年提出，对图像Image，记t为前景与背景的分割阈值，
前景点数占图像比例为w0，平均灰度为u0；背景点数占图像比例为w1，平均灰度为u1。
图像的总平均灰度为：u=w0*u0+w1*u1。从最小灰度值到最大灰度值遍历t，
当t使得值g=w0*(u0-u)2+w1* (u1-u)2 最大时t即为分割的最佳阈值。
 直接应用大津法计算量较大，因此我们在实现时采用了等价的公式g=w0*w1*(u0-u1)2。
部分计算过程如下'''
def OtsuGrayThreshold(image):
# 创建Hist
    #image = Image.open("4.f")
    image.convert('L')
    hist = image.histogram()
print  len(hist)
print hist

# 开始计算
    # 计算总亮度
    totalH = 0
for h in range(0,256):
        v =hist[h]
if v == 0 : continue
        totalH += v*h
print h,totalH,v*h


    width  = image.size[0]
    height = image.size[1]
    total  = width*height
print width
print height

#print "总像素:%d;总亮度:%d平均亮度:%0.2f"%(total,totalH,totalH/total)
    
    v = 0

    gMax = 0.0
    tIndex = 0

# temp
    n0Acc = 0.0
    n1Acc = 0.0
    n0H   = 0.0
    n1H   = 0.0
for t in range(1,255):
        v = hist[t-1]
if v == 0: continue

        n0Acc += v          #灰度小于t的像素的数目
        n1Acc = total - n0Acc #灰度大于等于t的像素的数目
        n0H += (t-1)*v          #灰度小于t的像素的总亮度
        n1H = totalH - n0H  #灰度大于等于t的像素的总亮度

if n0Acc > 0 and n1Acc > 0:
            u0 = n0H/n0Acc # 灰阶小于t的平均灰度
            u1 = n1H/n1Acc # 灰阶大于等于t的平均灰度
            w0 = n0Acc/total # 灰阶小于t的像素比例
            w1 = 1.0-w0      # 灰阶大于等于t的像素的比例
            uD = u0-u1
            g = w0 * w1 * uD * uD
#            print 'g=',g
            #if debug > 2: print "t=%3d; u0=%.2f,u1=%.2f,%.2f;n0H=%d,n1H=%d; g=%.2f"\
            #   %(t,u0,u1,u0*w0+u1*w1,n0H,n1H,g)
#            print t,u0,u1,w0,w1,g
            if gMax < g:
                gMax   = g
                tIndex = t

#    if debug >0 : print "gMaxValue=%.2f; t = %d ; t_inv = %d"\
#               %(gMax,tIndex,255-tIndex)            
#    print tIndex 
    return tIndex

def OtsuGray(image):
    otsuthreshold=OtsuGrayThreshold(image)
    bim = fxxBinary.fxxBinary(image, otsuthreshold)
return bim


if __name__ =='__main__':
    image = Image.open('roi.jpg')
    image = image.filter(ImageFilter.EDGE_ENHANCE())
    image = image.filter(ImageFilter.MedianFilter())
    image.show()
#enhancer = ImageEnhance.Contrast(im)
    otsu=OtsuGray(image)
    otsu.show()

4.边缘提取

#sobel算子 #perwit算子 #isotropic算子

# -*- coding: utf-8 -*-
from PIL import Image

def SobelSharp(image):  #sobel算子
    SobelX = [-1,0,1,-2,0,2,-1,0,1]
    SobelY = [-1,-2,-1,0,0,0,1,2,1]
    iSobelSharp = sharp(image,SobelX,SobelY)
return iSobelSharp

def PrewittSharp(image):#perwit算子
    PrewittX = [1,0,-1,1,0,-1,1,0,-1]
    PrewittY = [-1,-1,-1,0,0,0,1,1,1]
    iPrewittSharp = sharp(image,PrewittX,PrewittY)
return iPrewittSharp

def IsotropicSharp(image):#isotropic算子
    IsotropicX = [1,0,-1,1.414,0,-1.414,1,0,-1]
    IsotropicY = [-1,-1.414,-1,0,0,0,1,1.414,1]
    iIsotropicSharp = sharp(image,IsotropicX,IsotropicY)
return iIsotropicSharp

def sharp(image,arrayX,arrayY):
    w = image.size[0]
    h = image.size[1]
    size = (w,h)

    iSharpImage = Image.new('L', size)
    iSharp = iSharpImage.load()
    image = image.load()

    tmpX = [0]*9
    tmpY = [0]*9
for i in range(1,h-1):
for j in range(1,w-1):
for k in range(3):
for l in range(3):
                    tmpX[k*3+l] = image[j-1+l, i-1+k]*arrayX[k*3+l]
                    tmpY[k*3+l] = image[j-1+l, i-1+k]*arrayY[k*3+l]
            iSharp[j,i] = abs(sum(tmpX))+abs(sum(tmpY))
return iSharpImage


if __name__ =='__main__':
    image = Image.open('roi_open.jpg')
    iSobelSharp = SobelSharp(image)
    iSobelSharp.show()
    iPrewittSharp = PrewittSharp(image)
    iPrewittSharp.show()
    iIsotropicSharp = IsotropicSharp(image)
    iIsotropicSharp.show()

5.形态学开闭运算

   开闭运算

# -*- coding: utf-8 -*-
'''
Morphplogical operations such as expansion and erosion

'''
import operator, Image, ImageFilter

def  Corrode(image):
    w = image.size[0]
    h = image.size[1]
    size = (w, h)

    iCorrodeImage = Image.new('L',size)
    iCorrode = iCorrodeImage.load()
    image = image.load()


    kH = range(5)+range(h-5,h)
    kW = range(5)+range(w-5,w)
for i in range(h):
for j in range(w):
if i in kH or j in kW:
                iCorrode[j,i] = 0
elif image[j,i] == 255:
                iCorrode[j,i] = 255
else:
                a = []
for k in range(10):
for l in range(10):
                        a.append(image[j-5+l, i-5+k])
if max(a) == 255:
                    iCorrode[j, i] = 255
else:
                    iCorrode[j,i] = 0
return iCorrodeImage

def Expand(image):
    w = image.size[0]
    h = image.size[1]
    size = (w,h)

    iExpandImage = Image.new('L',size)
    iExpand = iExpandImage.load()
    image = image.load()

for i in range(w):
for j in range(h):
            iExpand[i,j] = 255
for i in range(h):
for j in range(w):
if image[j,i] == 0:
for k in range(10):
for l in range(10):
if -1<(i-5+k)<h and -1<(j-5+l)<w:
                            iExpand[j-5+l, i-5+k] = 0
return iExpandImage         

#测试代码
if __name__ == '__main__':
    image = Image.open('fc.bmp')
    iCorrodeImage = Corrode(image)
    iCorrodeImage.show()
    iExpandImage = Expand(image)
    iExpandImage.show()

6.骨架抽取细化

from PIL import Image

def Two(image):
    w = image.size[0]
    h = image.size[1]
    size = (w,h)
print size
    iTwoImage = Image.new('L', size)
    iTwo = iTwoImage.load()
    image = image.load()

for i in range(h):
for j in range(w):
#print image[j, i]
            iTwo[j,i] = 0 if image[j,i] > 200 else 255
return iTwoImage


def VThin(image,array):
    h = image.size[1]
    w = image.size[0]
    image = image.load()

    NEXT = 1
for i in range(h):
for j in range(w):
if NEXT == 0:
                NEXT = 1
else:
                M = image[j-1,i]+image[j,i]+image[j+1,i] if 0<j<w-1 else 1
if image[j, i] == 0  and M != 0:                  
                    a = [0]*9
for k in range(3):
for l in range(3):
if -1<(i-1+k)<h and -1<(j-1+l)<w and image[j-1+l, i-1+k]==255:
                                a[k*3+l] = 1
                    sum = a[0]*1+a[1]*2+a[2]*4+a[3]*8+a[5]*16+a[6]*32+a[7]*64+a[8]*128
                    image[j,i] = array[sum]*255
if array[sum] == 1:
                        NEXT = 0
return image

def HThin(image,array):
    h = image.size[1]
    w = image.size[0]
    image = image.load()

    NEXT = 1
for j in range(w):
for i in range(h):
if NEXT == 0:
                NEXT = 1
else:
                M = image[j, i-1]+image[j,i]+image[j, i+1] if 0<i<h-1 else 1   
if image[j, i] == 0 and M != 0:                  
                    a = [0]*9
for k in range(3):
for l in range(3):
if -1<(i-1+k)<h and -1<(j-1+l)<w and image[j-1+l, i-1+k]==255:
                                a[k*3+l] = 1
                    sum = a[0]*1+a[1]*2+a[2]*4+a[3]*8+a[5]*16+a[6]*32+a[7]*64+a[8]*128
                    image[j,i] = array[sum]*255
if array[sum] == 1:
                        NEXT = 0
return image

def Thining(image,array, num=10):
    iThinnedImage = image.copy()
for i in range(num):
        VThin(iThinnedImage,array)
        HThin(iThinnedImage,array)
return iThinnedImage

def Thin(image): 
    array = [0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\
1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1,\
             0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\
1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1,\
1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,\
             0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
1,1,0,0,1,1,0,0,1,1,0,1,1,1,0,1,\
             0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,\
             0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\
1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,1,\
             0,0,1,1,0,0,1,1,1,1,0,1,1,1,0,1,\
1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,\
1,1,0,0,1,1,0,0,0,0,0,0,0,0,0,0,\
1,1,0,0,1,1,1,1,0,0,0,0,0,0,0,0,\
1,1,0,0,1,1,0,0,1,1,0,1,1,1,0,0,\
1,1,0,0,1,1,1,0,1,1,0,0,1,0,0,0]

    image = image.convert('L')
#image.show()
    iTwo = Two(image)
#iTwo.show()
    iThin = Thining(iTwo, array)
#iThin.show()
    iRes = iThin.point(lambda i: 255-i)
return iRes



if __name__ == '__main__':
    image = Image.open('roi_edge.bmp')
    iThinImage = Thin(image)
    iThinImage.save('roi_pre.bmp')
    iThinImage.show()