# 打開一個新文件，將其命名為 watershed.py ，然后插入以下代碼：

# 導入必要的包
from skimage.feature import peak_local_max
from skimage.morphology import watershed
from scipy import ndimage
import numpy as np
import argparse
import imutils
import cv2

# 構(gòu)造參數(shù)解析并解析參數(shù)
ap = argparse.ArgumentParser()
# ap.add_argument("-i", "--image", default="HFOUG.jpg", help="path to input image")
ap.add_argument("-i", "--image", default="watershed_coins_01.jpg", help="path to input image")
args = vars(ap.parse_args())
# 加載圖像并執(zhí)行金字塔均值偏移濾波以輔助閾值化步驟
image = cv2.imread(args["image"])
shifted = cv2.pyrMeanShiftFiltering(image, 21, 51)
cv2.imshow("Input", image)
# 將圖像轉(zhuǎn)換為灰度，然后應(yīng)用大津閾值
gray = cv2.cvtColor(shifted, cv2.COLOR_BGR2GRAY)
thresh = cv2.threshold(gray, 0, 255, cv2.THRESH_BINARY | cv2.THRESH_OTSU)[1]
cv2.imshow("Thresh", thresh)

# 計算從每個二進制圖像中的像素到最近的零像素的精確歐氏距離，然后找出這個距離圖中的峰值
D = ndimage.distance_transform_edt(thresh)

# 可視化距離函數(shù)
D_show = cv2.normalize(D, None, 0, 1, cv2.NORM_MINMAX)
# print(np.max(D_show))
cv2.imshow("D_show", D_show)


# 以坐標列表(indices=True)或布爾掩碼(indices=False)的形式查找圖像中的峰值。峰值是2 * min_distance + 1區(qū)域內(nèi)的局部最大值。
# (即峰值之間至少相隔min_distance)。此處我們將確保峰值之間至少有20像素的距離。
localMax = peak_local_max(D, indices=False, min_distance=20, labels=thresh)
# 可視化localMax
temp = localMax.astype(np.uint8)
cv2.imshow("localMax", temp * 255)
# 使用8-連通性對局部峰值進行連接成分分析，然后應(yīng)用分水嶺算法
# scipy.ndimage.label(input, structure=None, output=None)
# input ：待標記的數(shù)組對象。輸入中的任何非零值都被視為待標記對象，零值被視為背景。
# structure：定義要素連接的結(jié)構(gòu)化元素。對于二維數(shù)組。默認是四連通， 此處選擇8連通
#
markers = ndimage.label(localMax, structure=np.ones((3, 3)))[0]  # [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

# 可視化markers
temp_markers = markers.astype(np.uint8)
cv2.imshow("temp_markers", temp_markers * 20)

# 由于分水嶺算法假設(shè)我們的標記代表距離圖中的局部最小值（即山谷），因此我們?nèi)?D 的負值。
labels = watershed(-D, markers, mask=thresh)
print("[INFO] {} unique segments found".format(len(np.unique(labels)) - 1))

# 循環(huán)遍歷分水嶺算法返回的標簽
for label in np.unique(labels):
    # 0表示背景，忽略它
    if label == 0:
        continue
    # 否則，為標簽區(qū)域分配內(nèi)存并將其繪制在掩碼上
    mask = np.zeros(gray.shape, dtype="uint8")
    mask[labels == label] = 255
    # 在mask中檢測輪廓并抓取最大的一個
    cnts = cv2.findContours(mask.copy(), cv2.RETR_EXTERNAL,
                            cv2.CHAIN_APPROX_SIMPLE)
    cnts = imutils.grab_contours(cnts)
    c = max(cnts, key=cv2.contourArea)
    # 在物體周圍畫一個圓
    ((x, y), r) = cv2.minEnclosingCircle(c)
    cv2.circle(image, (int(x), int(y)), int(r), (0, 255, 0), 2)
    cv2.putText(image, "#{}".format(label), (int(x) - 10, int(y)),
                cv2.FONT_HERSHEY_SIMPLEX, 0.6, (0, 0, 255), 2)
# 顯示輸出圖像
cv2.imshow("Output", image)
cv2.waitKey(0)

2.Watershed與random walker分割對比

示例比較了兩種分割方法，以分離兩個相連的磁盤：分水嶺算法和隨機游走算法。兩種分割方法都需要種子，即明確屬于某個區(qū)域的像素。在這里，到背景的距離圖的局部最大值被用作種子。

import numpy as np
from skimage.segmentation import watershed
from skimage.feature import peak_local_max
from skimage import measure
from skimage.segmentation import random_walker
import matplotlib.pyplot as plt
from scipy import ndimage
# import cv2

# Generate an initial image with two overlapping circles
x, y = np.indices((80, 80))
x1, y1, x2, y2 = 28, 28, 44, 52
r1, r2 = 16, 20
mask_circle1 = (x - x1) ** 2 + (y - y1) ** 2 < r1 ** 2
mask_circle2 = (x - x2) ** 2 + (y - y2) ** 2 < r2 ** 2
image = np.logical_or(mask_circle1, mask_circle2)
# Now we want to separate the two objects in image
# Generate the markers as local maxima of the distance
# to the background
distance = ndimage.distance_transform_edt(image)

D_show = distance/np.max(distance)

# D_show = cv2.normalize(distance, None, 0, 1, cv2.NORM_MINMAX)

local_maxi = peak_local_max(
    distance, indices=False, footprint=np.ones((3, 3)), labels=image)
markers = measure.label(local_maxi)
labels_ws = watershed(-distance, markers, mask=image)

markers[~image] = -1
labels_rw = random_walker(image, markers)

plt.figure(figsize=(12, 3.5))
plt.subplot(141)
plt.imshow(image, cmap='gray', interpolation='nearest')
plt.axis('off')
plt.title('image')
plt.subplot(142)
plt.imshow(D_show, cmap='Spectral',interpolation='nearest')
plt.axis('off')
plt.title('distance map')
plt.subplot(143)
plt.imshow(labels_ws, cmap='Spectral', interpolation='nearest')
plt.axis('off')
plt.title('watershed segmentation')
plt.subplot(144)
plt.imshow(labels_rw, cmap='Spectral', interpolation='nearest')
plt.axis('off')
plt.title('random walker segmentation')

plt.tight_layout()
plt.show()