import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
def np_from_img(fname):
img = Image.open(fname)
print img.format, img.size, img.mode
img = img = img.convert('I')
#print img.format, img.size, img.mode
return np.asarray(img)
#########
img = np_from_img('E:/tmp/python/img/2.bmp')
class ImageSegmentationMcmPde(object):
def __init__(self):
self.iter = 0
self.curv_index = None
self.curve_img = None
self.num = 0
self.init_edge = None
self.m = 0
self.n = 0
self.U = None
def init_config(self,img):
self.m = int(img.shape[0])
self.n = int(img.shape[1])
self.U = np.copy(img)
print type(self.m)
self._cal_border(True)
self._cal_dist(True)
def _cal_border(self,init=False):
m = self.m
n = self.n
border_min = 0
border_max = 0
if init == True:
border_min = 5
border_max = 120
curve_img = np.zeros((m,n))
num = 0
curv_index = []
for i in range(1,m-1):
for j in range(1,n-1):
if self.U[i][j] < border_min and ( self.U[i - 1][j] > border_max or self.U[i + 1][j] > border_max or self.U[i][j -1] > border_max or self.U[i][j + 1] > border_max):
is_exist = False
for k in range(0,num):
if curv_index[k][0] == i and curv_index[k][1] == j:
is_exist = True
break
else:
pass
if is_exist == False:
curv_index.append((i,j))
curve_img[i][j] = 255
num = num + 1
self.curv_index = curv_index
self.curve_img = curve_img
self.num = num
if init == True:
self.init_edge = np.copy(curve_img);
def _cal_dist(self, init=False):
m = self.m
n = self.n
num = self.num
bmin = 0
if init == True:
bmin = 5
curv_index = self.curv_index
U = np.zeros((m,n))
for i in range(m):
print 'iter(%s,%s)' % (self.iter,i)
for j in range(n):
mmin = 100000
for k in range(num):
sq = (i - curv_index[k][0]) * (i - curv_index[k][0]) + (j - curv_index[k][1]) * (j - curv_index[k][1])
sq = np.sqrt(sq)
if sq < mmin:
mmin = sq
U[i][j] = mmin
if self.U[i][j] < bmin:
U[i][j] *= -1
self.U = U
def do_iter(self):
dt = 0.1
U = self.U
m = self.m
n = self.n
# display_img = init_curve_img
diffx1 = np.concatenate((np.arange(1,n),[n-1]))
diffx2 = np.concatenate(([0],np.arange(0,n-1)))
diffy1 = np.concatenate((np.arange(1,m),[m-1]))
diffy2 = np.concatenate(([0],np.arange(0,m-1)))
#U_x = (U[:,[2:nny ,nny]]-U[:,[1 1:nny-1]]/2
U_x = (U[:, diffx1] - U[:, diffx2]) /2
#U_y = (U[[2:nnx nnx],:]-U[[1 1:nnx-1],:]]/2
U_y = (U[diffy1,:] - U[diffy2,:])/2
# U_xx = U[:,[2:nny nny]]+U[:,[1 1:nny-1]]-2*U
# U_yy = U[[2:nnx nnx],:]+U[[1 1:nnx-1],:]-2*U
U_xx = U[:,diffx1] + U[:,diffx2] - 2 * U
U_yy = U[diffy1,:] + U[diffy2,:] - 2 * U
# Dp = U([2:nnx nnx],[2:nny nny])+U([1 1:nnx-1],[1 1:nny-1])
# Dm = U([1 1:nnx-1],[2:nny nny])+U([2:nnx nnx],[1 1:nny-1])
#Dp = U[diffy2,diffx2] + U[diffy2,diffx2]
Dp = U[diffy1,:][:,diffx1] + U[diffy2,:][:,diffx2]
Dm = U[diffy2,:][:,diffx1] + U[diffy1,:][:,diffx2]
U_xy = (Dp-Dm)/4
# Num = U_xx.*(U_y.^2)-2*U_x.*U_y.*U_xy+U_yy.*(U_x.^2)
# Den = U_x.^2+U_y.^2
# I_t = Num./(Den+eps)
Num = U_xx * U_y * U_y - 2 * U_x * U_y * U_xy + U_yy * U_x * U_x
Den = U_x * U_x + U_y * U_y
I_t = Num / (Den + 0.000001)
U = U + dt * I_t
self.iter = self.iter + 1
if self.iter % 1000 == 0:
self._cal_border()
self._cal_dist()
self.init_edge = self.init_edge + self.curve_img
#plt.matshow(self.init_edge)
#plt.show()
def do_test(self):
for i in range(5000):
self.do_iter()
plt.matshow(self.init_edge)
plt.show()
###################
if __name__ == '__main__':
imgPde = ImageSegmentationMcmPde()
imgPde.init_config(img)
imgPde.do_test()
参考
<<图像处理的偏微分方法>>
<<Geometric partial differential equations and image analysis (Sapiro)>>