前言
形态学建筑物指数MBI通过建立建筑物的隐式特征和形态学算子之间的关系进行建筑物的提取[1]。
原理
上图源自[2]。
实验数据
简单找了一张小图片：
test.jpg
代码
为了支持遥感图像，读写数据函数都是利用GDAL写的。
import numpy as np
import gdal# 读取tif数据集
def readTif(fileName, xoff = 0, yoff = 0, data_width = 0, data_height = 0):dataset = gdal.Open(fileName)if dataset == None:print(fileName + "文件无法打开")# 栅格矩阵的列数width = dataset.RasterXSize # 栅格矩阵的行数height = dataset.RasterYSize # 波段数bands = dataset.RasterCount # 获取数据if(data_width == 0 and data_height == 0):data_width = widthdata_height = heightdata = dataset.ReadAsArray(xoff, yoff, data_width, data_height)# 获取仿射矩阵信息geotrans = dataset.GetGeoTransform()# 获取投影信息proj = dataset.GetProjection()return width, height, bands, data, geotrans, proj# 保存tif文件函数
def writeTiff(im_data, im_geotrans, im_proj, path):if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelif len(im_data.shape) == 2:im_data = np.array([im_data])im_bands, im_height, im_width = im_data.shape#创建文件driver = gdal.GetDriverByName("GTiff")dataset = driver.Create(path, int(im_width), int(im_height), int(im_bands), datatype)if(dataset!= None):dataset.SetGeoTransform(im_geotrans) #写入仿射变换参数dataset.SetProjection(im_proj) #写入投影for i in range(im_bands):dataset.GetRasterBand(i+1).WriteArray(im_data[i])del dataset接下来就是就算MBI，代码注释很详细，也可以对着原理来看。
from skimage.morphology import square, white_tophat
from skimage.transform import rotate# 计算MBI
# s_min: 结构元素大小最小值
# s_max: 结构元素大小最大值
# delta_s: 颗粒测定的间隔
def CalculationMBI(filePath, MBIPath, s_min, s_max, delta_s):# 读取图像的相关信息width, height, bands, image, geotrans, proj = readTif(filePath)# 多光谱带的最大值对应于具有高反射率的特征->取光谱带最大值作为后续计算数据gray = np.max(image, 0)# 为消除白帽边缘效应，进行边缘补零gray = np.pad(gray, ((s_min, s_min), (s_min, s_min)), 'constant', constant_values=(0, 0))# 形态学剖面集合MP_MBI_list = []# 差分形态学剖面DMP集合DMP_MBI_list = []# 计算形态学剖面for i in range(s_min, s_max + 1, 2 * delta_s):print("s = ", i)# 大小为i×i的单位矩阵SE_intermediate = square(i)# 只保留中间一行为1,其他设置为0SE_intermediate[ : int((i - 1) / 2), :] = 0SE_intermediate[int(((i - 1) / 2) + 1) : , :] = 0# SE_intermediate表示结构元素，用于设定局部区域的形状和大小# 旋转0 45 90 135°for angle in range(0, 180, 45):SE_intermediate = rotate(SE_intermediate, angle, order = 0, preserve_range = True).astype('uint8')# 多角度形态学白帽重构MP_MBI = white_tophat(gray, selem = SE_intermediate)MP_MBI_list.append(MP_MBI)# 计算差分形态学剖面DMPfor j in range(4, len(MP_MBI_list), 1):# 差的绝对值DMP_MBI = np.absolute(MP_MBI_list[j] - MP_MBI_list[j - 4])DMP_MBI_list.append(DMP_MBI)# 计算MBIMBI = np.sum(DMP_MBI_list, axis = 0) / (4 * (((s_max - s_min) / delta_s) + 1))# 去除多余边缘结果MBI = MBI[s_min : MBI.shape[0] - s_min, s_min : MBI.shape[1] - s_min]# 写入文件writeTiff(MBI, geotrans, proj, MBIPath)# 原图像
filePath = r"test.jpg"
# MBI结果
MBIPath = r"test_mbi.jpg"
# 建筑物提取结果
buildingPath = r"test_building.jpg"
# 结构元素大小最小值
s_min = 3
# 结构元素大小最大值
s_max = 20
# 测定的间隔
delta_s = 1
# 计算MBI
CalculationMBI(filePath, MBIPath, s_min, s_max, delta_s)
test_mbi.jpg
MBI计算出来了以后，我们就要取阈值来提取建筑物了，阈值可以手动设置，也可以用算法自动求出阈值，这里我们采用OTSU算法[3]。
from skimage.filters import threshold_otsudef BuildingExtraction_otsu(MBIPath, buildingPath):width, height, bands, image, geotrans, proj = readTif(MBIPath)thresh = threshold_otsu(image) #返回一个阈值image[image>thresh] = 255image[image<=thresh] = 0image = image.astype(np.uint8)writeTiff(image, geotrans, proj, buildingPath)# otsu自动计算阈值提取建筑物
BuildingExtraction_otsu(MBIPath, buildingPath)
test_building.jpg
目视对照一下的话，感觉效果还不错。
参考
-
^Huang X and Zhang L. 2011. A multidirectional and multiscale morphological index for automatic building extraction from multispectral geoeye-1 imagery. Photogrammetric Engineering and Remote Sensing, 77(7), 721-732. [DOI: 10.14358/PERS.77.7.721]
-
^魏旭,高小明,岳庆兴,郭正胜.一种结合MBI和SLIC算法的遥感影像建筑物提取方法[J].测绘与空间地理信息,2019,42(10):100-103.
-
^otsu（大津算法）-百度百科 https://baike.baidu.com/item/otsu/16252828?fr=aladdin
来源：应用推广部
供稿：技术研发部
编辑：方梅
