#!/usr/bin/env python
# -*- coding: utf-8 -*-
# @Time : 2022/4/16 9:53
# @Author : suxiaorui
import matplotlib.pyplot as plt
import cv2
import numpy as np
def sdf2circle(row,col,ic,jc,r):
'''
:param row: 行数
:param col: 列数
:param ic: 圆心
:param jc: 圆心
:param r: 圆的半径
:return: 到圆的距离
'''
X = [i for i in range(col)]
Y = [i for i in range(row)]
X,Y = np.meshgrid(X,Y)
f = np.sqrt((X-jc)**2 + (Y-ic)**2) - r
return f
def plotLevelSet(u,zLevel, style):
return plt.contour(u, levels=[zLevel], colors=[style])
def BoundMirrorExpand(A):
m,n = A.shape[0],A.shape[1]
B = np.zeros((m+2,n+2))
B[1:m+1,1:n+1] = A
# 处理边界
for i in range(n+1):
B[0,i] = B[2,i]
B[m+1,i] = B[m-1,i]
for i in range(m+2):
B[i,0] = B[i,2]
B[i,n+1] = B[i,n-1]
return B
def BoundMirrorEnsure(A):
m,n = A.shape[0],A.shape[1]
if m < 3 or n < 3:
print("either the number of rows or columns is smaller than 3")
B = A
for i in range(n):
B[0, i] = B[2, i]
B[m - 1, i] = B[m - 3, i]
for i in range(m):
B[i, 0] = B[i, 2]
B[i, n - 1] = B[i, n - 3]
return B
def Delta(phi, epsilon):
return epsilon/(np.pi*(phi**2+epsilon**2))
def Heaviside(phi,epsilon):
return 0.5 * (1 + (2 / np.pi) * np.arctan(phi/epsilon))
def curvature(f):
'''
计算 % K=div(Df/|Df|) = (fxx*fy^2+fyy*fx^2-2*fx*fy*fxy)/(fx^2+fy^2)^(3/2)
'''
# 这里使用的是中心差分来计算的
Ix, Iy = np.gradient(f)
normDu = np.sqrt(Ix**2 + Iy**2 + 1e-10)
Nx = Ix / normDu
Ny = Iy / normDu
[nxx, junk] = np.gradient(Nx)
[junk, nyy] = np.gradient(Ny)
K = nxx + nyy
return K
def binaryfit(phi,U,epsilon):
'''
:param phi: 水平集函数
:param U: 输入的图像
:param epsilon: 参数
:return: 两个值 C1 在区域内部 C2在区域外部
'''
H = Heaviside(phi, epsilon)
a = H * U
numer_1 = np.sum(a)
denom_1 = np.sum(H)
C1 = numer_1 / denom_1
b = (1 - H) * U
numer_2 = np.sum(b)
c = 1 - H
denom_2 = np.sum(c)
C2 = numer_2 / denom_2
return [C1,C2]
def BoundMirrorShrink(A):
m,n = A.shape[0],A.shape[1]
B = A[1:m-1,1:n-1]
return B
def evolution_cv(I, phi0, mu, nu, lambda_1, lambda_2, delta_t, epsilon, numIter):
'''
:param I: 输入的图像
:param phi0: 要更新的水平集函数
:param mu: 长期项的权重
:param nu: 面积项的权重
:param lambda_1: 拟合项权重1
:param lambda_2: 拟合项权重1
:param delta_t: 时间步长
:param epsilon: 用来计算Heaviside 和 dirac 函数的参数
:param numIter: 迭代次数
:return: 更细的水平集函数
'''
I = BoundMirrorExpand(I)
phi = BoundMirrorExpand(phi0)
for k in range(numIter):
phi = BoundMirrorEnsure(phi)
delta_h = Delta(phi, epsilon)
Curv = curvature(phi)
[C1, C2] = binaryfit(phi, I, epsilon)
# 更新phi
phi = phi + delta_t * delta_h * (mu*Curv-nu-lambda_1*(I-C1)**2+lambda_2*(I-C2)**2)
phi = BoundMirrorShrink(phi)
return phi
Image = cv2.imread('1.bmp')
Image = cv2.cvtColor(Image,cv2.COLOR_BGR2GRAY)
U = np.array(Image)
row,col = U.shape[0],U.shape[1]
ic = row / 2
jc = col / 2
r = 30
phi_0 = sdf2circle(row,col,ic,jc,r)
delta_t = 0.1
lambda_1 = 1
lambda_2 = 1
nu = 0
h = 1
epsilon = 1
mu = 0.03 * 255 * 255
I = np.array(U,dtype=np.float64)
phi = phi_0
img = cv2.cvtColor(Image,cv2.COLOR_BGR2RGB)
plt.figure(1)
plt.imshow(img)
# plotLevelSet(phi,0,'r')
# plt.show()
plt.draw()
numIter = 1
for k in range(0,1001):
phi = evolution_cv(I,phi,mu,nu,lambda_1,lambda_2,delta_t,epsilon,numIter)
if k % 200 == 0:
print('第'+str(k)+'次：')
plt.imshow(img)
plotLevelSet(phi,0,'r')
plt.draw()
plt.show()
plt.pause(0.01)
- 1
- 2
- 3
- 4
前往页