手背血管预处理--利用PIL进行数字图像的综合处理_
终于把这个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: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')
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()
posted on 2011-05-20 17:16 追求卓越 挑战极限 阅读(2136) 评论(4) 编辑 收藏 举报
