1. 定義
分水嶺算法(watershed algorithm)可以將圖像中的邊緣轉化為“山脈”，將均勻區域轉化為“山谷”，在這方面有助於分割目標。
分水嶺算法：是一種基於拓撲理論的數學形態學的分割方法。把圖像看作是測地學上的拓撲地貌，圖像中的每一個點像素值的灰度值表示該點的海拔高度，每一個局部極小值及其影響的區域稱為“集水盆”，集水盆的邊界可以看成分水嶺。在每一個局部極小值表面刺穿一個小孔，然后把整個模型慢慢浸入水中，隨着浸入的加深，每一個局部極小值的影響域慢慢的向外擴展，在兩個集水盆匯合處構建大壩，形成分水嶺。
迭代標注過程：
- 排序過程：對每個像素的灰度級進行從低到高的排序
- 淹沒過程：對每一個局部最小值在h階高度的影響域采用先進先出結構判斷及標注。
2.實現算法:watershed()函數
這些標記的值可以使用findContours（）函數和drawContours（）函數由二進制的掩模檢索出來
3.程序代碼：
import numpy as np import cv2 from osgeo import gdal, gdal_array import shapefile try: import Image import ImageDraw except: from PIL import Image, ImageDraw def tif_jpg(rasterfile): in_ds = gdal.Open(rasterfile) # 打開樣本文件 xsize = in_ds.RasterXSize # 獲取行列數 ysize = in_ds.RasterYSize bands = in_ds.RasterCount block_data = in_ds.ReadAsArray(0, 0, xsize, ysize).astype(np.float32) B = block_data[0, :, :] G = block_data[ 1,:, :] R = block_data[2,:, :] R1 = (R/np.max(R)*255).astype(np.int16) G1 = (G / np.max(G) * 255).astype(np.int16) B1 = (B / np.max(B) * 255).astype(np.int16) data2 = cv2.merge([R1,G1,B1]) return data2 def watershed(path,out): print("分水嶺分割") in_ds = gdal.Open(path) # 打開樣本文件 xsize = in_ds.RasterXSize # 獲取行列數 ysize = in_ds.RasterYSize bands = in_ds.RasterCount geotransform = in_ds.GetGeoTransform() projection = in_ds.GetProjectionRef() #tif轉jpg 非255通道轉換為255通道 img=tif_jpg(path).astype(np.uint8) # 轉換為灰度圖片 gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY) # canny邊緣檢測 函數返回一副二值圖，其中包含檢測出的邊緣。 canny = cv2.Canny(gray_img, 80,120) # 尋找圖像輪廓 返回修改后的圖像 圖像的輪廓 以及它們的層次 # canny, contours, hierarchy = cv2.findContours(canny, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) contours, hierarchy = cv2.findContours(canny, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE) # 32位有符號整數類型， marks = np.zeros(img.shape[:2], np.int32) # findContours檢測到的輪廓 imageContours = np.zeros(img.shape[:2], np.uint8) # 輪廓顏色 compCount = 0 index = 0 # 繪制每一個輪廓 for index in range(len(contours)): # 對marks進行標記，對不同區域的輪廓使用不同的亮度繪制，相當於設置注水點，有多少個輪廓，就有多少個輪廓 # 圖像上不同線條的灰度值是不同的，底部略暗，越往上灰度越高 marks = cv2.drawContours(marks, contours, index, (index, index, index), 1, 8, hierarchy) # 繪制輪廓，亮度一樣 imageContours = cv2.drawContours(imageContours, contours, index, (255, 255, 255), 1, 8, hierarchy) # 查看 使用線性變換轉換輸入數組元素成8位無符號整型。 markerShows = cv2.convertScaleAbs(marks) # cv2.imshow('imageContours',imageContours) # 使用分水嶺算法 marks = cv2.watershed(img, marks) driver = gdal.GetDriverByName('GTiff') outfile_lake = out + "\\" + "watershed_cut.tif" out_dataset = driver.Create(outfile_lake, xsize, ysize, 1, gdal.GDT_Float32) out_band1 = out_dataset.GetRasterBand(1) out_band1.WriteArray(marks) out_dataset.SetGeoTransform(geotransform) # 寫入仿射變換 out_dataset.SetProjection(projection) return outfile_lake if __name__ == "__main__": path = r"D:\data\實驗數據\fenlei2.tif" out = r"D:\data\實驗結果\分割結果" watershed(path, out)