基于自選種子點的區域生長（python）

Introduction

當單一的分割算法的提取結果出現欠分割問題時，需要考慮加入影像後處理，進行圖像同質區的合并。本博文将簡單介紹下如何利用區域生長法來進行該項處理。如下圖所示，最右邊是欠分割掩膜，中間是同質區合并後的結果。

PS:區域生長法的算法設計思路：

區域生長算法的設計主要由以下三點：生長種子點的确定，區域生長的條件，區域生長停止的條件。種子點的個數根據具體的問題可以選擇一個或者多個，并且根據具體的問題不同可以采用完全自動确定或者人機互動确定。

區域生長的條件實際上就是根據像素灰階間的連續性而定義的一些相似性準則，而區域生長停止的條件定義了一個終止規則，基本上，在沒有像素滿足加入某個區域的條件的時候，區域生長就會停止。在算法裡面，定義一個變量，最大像素灰階值距離reg_maxdist.當待加入像素點的灰階值和已經分割好的區域所有像素點的平均灰階值的差的絕對值小于或等于reg_maxdist時，該像素點加入到已經分割到的區域。相反，則區域生長算法停止。

在種子店1的4鄰域連通像素中，即2、3、4、5點，像素點5的灰階值與種子點的灰階值最接近，是以像素點5被加入到分割區域中，并且像素點5會作為新的種子點執行後面的過程。在第二次循環過程中，由于待分析圖像中，即2、3、4、6、7、8，像素7的灰階值和已經分割的區域（由1和5組成）的灰階均值10.5最接近，是以像素點7被加入到分割區域中。最後一幅圖，示意了區域生長的方向（由淺入深）。

code

# -*- coding:utf-8 -*-
import cv2
import numpy as np
import gdalTools
from scipy import ndimage as ndi
from skimage.morphology import remove_small_holes, closing, square, opening, remove_small_objects, watershed

import matplotlib.pyplot as plt
####################################################################################


#######################################################################################
class Point(object):
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def getX(self):
        return self.x

    def getY(self):
        return self.y


connects = [Point(-1, -1), Point(0, -1), Point(1, -1), Point(1, 0), Point(1, 1), Point(0, 1), Point(-1, 1),
            Point(-1, 0)]


#####################################################################################
# 計算兩個點間的歐式距離
def get_dist(image, seed_location1, seed_location2):
    l1 = image[seed_location1.x, seed_location1.y]
    l2 = image[seed_location2.x, seed_location2.y]
    count = np.sqrt(np.sum(np.square(l1 - l2)))
    return count


########################################################################################
def reginGrow(image, mask):
    im_shape = image.shape
    height = im_shape[0]
    width = im_shape[1]

    markers = ndi.label(mask, output=np.uint32)[0]
    unis = np.unique(markers)
    # 擷取種子點
    seed_list = []

    for uni in unis:
        if uni == 0:
            continue
        pointsX, pointsY = np.where(markers == uni)
        num_point = len(pointsX) // 4
        for i in [0, num_point * 1, num_point * 2, num_point * 3]:
            pointX, pointY = pointsX[i], pointsY[i]
            seed_list.append(Point(pointX, pointY))

    # 标記，判斷種子是否已經生長
    img_mark = np.zeros([height, width])

    T = 7.5  # 門檻值
    class_k = 1  # 類别
    # 生長一個類
    while (len(seed_list) > 0):
        seed_tmp = seed_list[0]
        # 将以生長的點從一個類的種子點清單中删除
        seed_list.pop(0)

        img_mark[seed_tmp.x, seed_tmp.y] = class_k

        # 周遊8鄰域
        for i in range(8):
            tmpX = seed_tmp.x + connects[i].x
            tmpY = seed_tmp.y + connects[i].y

            if (tmpX < 0 or tmpY < 0 or tmpX >= height or tmpY >= width):
                continue
            dist = get_dist(image, seed_tmp, Point(tmpX, tmpY))
            # 在種子集合中滿足條件的點進行生長
            if (dist < T and img_mark[tmpX, tmpY] == 0):
                img_mark[tmpX, tmpY] = class_k
                seed_list.append(Point(tmpX, tmpY))
    img_mark = img_mark + mask
    img_mark = remove_small_holes(img_mark.astype(np.uint8), 100)
    return np.where(img_mark > 0, 1, 0)


if __name__ == '__main__':

    #import Image
    im_proj, im_geotrans, im_width, im_height, im = gdalTools.read_img('data/image.tif')
    im = im.transpose((1, 2, 0))
    image = im.copy()
    _, _, _, _, mask = gdalTools.read_img('data/seed.tif')
    img_mark = reginGrow(image, mask)