Python--level set （水平集）和 chan-vese模型

level set :https://www.zhihu.com/question/22608763?sort=created

https://blog.csdn.net/xiangyong58/article/details/11876019

chan-vese模型(公式推导）：https://blog.csdn.net/zhangchen1003/article/details/48930377

水平集（CV模型）代码：

1. import cv2
2. from pylab import*
3. Image = cv2.imread('02.jpg', 1) # 读入原图
4. image = cv2.cvtColor(Image, cv2.COLOR_BGR2GRAY)
5. img = np.array(image, dtype=np.float64) # 读入到np的array中，并转化浮点类型
6. # 初始水平集函数
7. IniLSF = np.ones((img.shape[0], img.shape[1]), img.dtype)
8. IniLSF[300:320, 300:320] = -1
9. IniLSF = -IniLSF
10. # 画初始轮廓
11. Image = cv2.cvtColor(Image, cv2.COLOR_BGR2RGB)
12. plt.figure(1), plt.imshow(Image), plt.xticks([]), plt.yticks([]) # to hide tick values on X and Y axis
13. plt.contour(IniLSF, [0], color='b', linewidth=2) # 画LSF=0处的等高线
14. plt.draw(), plt.show(block=False)
15. def mat_math(intput, str):
16. output = intput
17. for i in range(img.shape[0]):
18. for j in range(img.shape[1]):
19. if str == "atan":
20. output[i, j] = math.atan(intput[i, j])
21. if str == "sqrt":
22. output[i, j] = math.sqrt(intput[i, j])
23. return output
24. # CV函数
25. def CV(LSF, img, mu, nu, epison, step):
26. Drc = (epison / math.pi) / (epison*epison + LSF*LSF)
27. Hea = 0.5*(1 + (2 / math.pi)*mat_math(LSF/epison, "atan"))
28. Iy, Ix = np.gradient(LSF)
29. s = mat_math(Ix*Ix+Iy*Iy, "sqrt")
30. Nx = Ix / (s+0.000001)
31. Ny = Iy / (s+0.000001)
32. Mxx, Nxx = np.gradient(Nx)
33. Nyy, Myy = np.gradient(Ny)
34. cur = Nxx + Nyy
35. Length = nu*Drc*cur
36. Lap = cv2.Laplacian(LSF, -1)
37. Penalty = mu*(Lap - cur)
38. s1 = Hea*img
39. s2 = (1-Hea)*img
40. s3 = 1-Hea
41. C1 = s1.sum() / Hea.sum()
42. C2 = s2.sum() / s3.sum()
43. CVterm = Drc*(-1 * (img - C1)*(img - C1) + 1 * (img - C2)*(img - C2))
44. LSF = LSF + step*(Length + Penalty + CVterm)
45. # plt.imshow(s, cmap ='gray'),plt.show()
46. return LSF
47. # 模型参数
48. mu = 1
49. nu = 0.003 * 255 * 255
50. num = 20
51. epison = 1
52. step = 0.1
53. LSF = IniLSF
54. for i in range(1, num):
55. LSF = CV(LSF, img, mu, nu, epison, step) # 迭代
56. if i % 1 == 0: # 显示分割轮廓
57. plt.imshow(Image), plt.xticks([]), plt.yticks([])
58. plt.contour(LSF, [0], colors='r', linewidth=2)
59. plt.draw(), plt.show(block=False), plt.pause(0.01)

为什么上传图片这么麻烦。