手背血管预处理--利用PIL进行数字图像的综合处理_ 手背血管预处理--利用PIL进行数

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

文章图片

（此图如有版权问题请与我联系）
ROI：
【手背血管预处理--利用PIL进行数字图像的综合处理_】
文章图片

最后的细化结果:

文章图片

网上有很多人采取的是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:ChenxofHit@gmail.com 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 ' ) 3Otsu阈值分割
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=https://www.it610.com/article/%.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() 转载于:https://www.cnblogs.com/ChenxofHit/archive/2011/05/20/2052182.html