This algorithm was introduced before , But it works IDL Language , You can see ： Ice lake extraction algorithm based on bimodal threshold segmentation （IDL Language implementation ）, Recently, due to the need to re-use this algorithm , Just use it python Language implements the following , Let's go straight to the code .
PS, One thing to note here is that , Bimodal threshold segmentation can be understood as weighted average value related to variance , So the denominator should be variance , Not the average .
Add another article to extract ice lake from random forest ： Ice lake extraction algorithm based on random forest （python Language implementation ）
def bimodal(img, save_path, initial_threshold=0.1):
if not os.path.exists(save_path):
os.makedirs(save_path)
##################### Ice lake crude extraction #####################
# Calculate the global NDWI
ndwi = (img[:, :, 4] - img[:, :, 2])/(img[:, :, 4] + img[:, :, 2] - 0.00000001)
ndwi = np.clip(ndwi, -1.0, 1.0)
initial_mask = np.zeros_like(ndwi)
for i in range(ndwi.shape[0]):
for j in range(ndwi.shape[1]):
if ndwi[i, j] > initial_threshold and img[i, j, 4] < 12500 and img[i, j, 5] < 7500:
initial_mask[i, j] = 1
# Filter to remove noise , Less than 5 A region of pixels is removed directly
whole_result = morphology.remove_small_objects(initial_mask.astype('int64'), min_size=5, connectivity=2)
# Filling cavity
whole_result = ndimage.binary_fill_holes(whole_result, structure=None, output=None, origin=0)
##################### Ice lake extract #####################
# Statistics of connected regions
temp1 = measure.label(whole_result, connectivity=2)
num_area = temp1.max()
final_mask = np.zeros_like(initial_mask)
for i in range(1, num_area):
# Build buffers one by one
temp2 = (temp1 == i) # Find No i Regions
num_pixel_before = sum(sum(temp2)) # Count the number of pixels in this area
if num_pixel_before < 5:
continue
# Inflate build buffer , Until the area is the same 2 Double up
temp3 = temp2
num_pixel_after = num_pixel_before
while num_pixel_after < 2 * num_pixel_before:
temp3 = morphology.dilation(temp3, morphology.square(3))
num_pixel_after = sum(sum(temp3))
# If the area is too large , We must consider whether it is an ice lake
if num_pixel_after > 10000:
break
# Calculate the foreground and background in the buffer
lake_ndwi, land_ndwi = [], []
for i in range(temp2.shape[0]):
for j in range(temp2.shape[1]):
if temp3[i, j] == 1:
lake_ndwi.append(ndwi[i, j])
if temp2[i, j] == 0:
land_ndwi.append(ndwi[i, j])
u_lake = sum(lake_ndwi) / len(lake_ndwi)
sigma_lake = math.sqrt(sum((lake_ndwi - u_lake)**2) / len(lake_ndwi))
u_land = sum(land_ndwi) / len(land_ndwi)
sigma_land = math.sqrt(sum((land_ndwi - u_land)**2) / len(land_ndwi))
bimodal_threshold = (u_land * sigma_lake + u_lake * sigma_land) / (sigma_lake + sigma_land)
print('threshold = %f' % bimodal_threshold)
# Split the results in the buffer
temp4 = ((ndwi * temp3) > bimodal_threshold)
# preservation
for i in range(temp4.shape[0]):
for j in range(temp4.shape[1]):
if temp4[i, j] == 1:
final_mask[i, j] = 1
io.imsave('threshold.png', final_mask)Direct reading Landsat-8 Remote sensing image , And then call this function , The result of bimodal threshold segmentation can be generated . Here I also assisted NIR<15000 and SWIR<7500 For further screening , This will filter out the effects of glaciers as much as possible .
Directly above ：
You can see , In fact, the bimodal threshold segmentation algorithm , If not assisted DEM Under the circumstances , In fact, the extraction effect of ice lake is not very good , At the edge of the ice lake , And glaciers 、 Rivers will be extracted . As mentioned before , In fact, the ice lake and background sometimes do not conform to the characteristics of bimodal distribution , So in the process of mapping large-scale ice lakes , The accuracy of this method is not high .
