MCM(平均曲率演化)

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)>>