119 lines
3.8 KiB
Python
119 lines
3.8 KiB
Python
from util import *
|
||
from osgeo import gdal
|
||
import argparse
|
||
|
||
|
||
@timeit
|
||
def otsu(img, max_value, data_water_mask, ignore_value=0, foreground=1, background=0):
|
||
height = img.shape[0]
|
||
width = img.shape[1]
|
||
|
||
hist = np.zeros([max_value], np.float32)
|
||
|
||
# 计算直方图
|
||
invalid_counter = 0
|
||
for i in range(height):
|
||
for j in range(width):
|
||
if img[i, j] == ignore_value or img[i, j] < 0 or data_water_mask[i, j] == 0:
|
||
invalid_counter = invalid_counter + 1
|
||
continue
|
||
|
||
hist[img[i, j]] += 1
|
||
|
||
hist /= (height * width - invalid_counter)
|
||
|
||
threshold = 0
|
||
deltaMax = 0
|
||
# 遍历像素值，计算最大类间方差
|
||
for i in range(max_value):
|
||
wA = 0
|
||
wB = 0
|
||
uAtmp = 0
|
||
uBtmp = 0
|
||
uA = 0
|
||
uB = 0
|
||
u = 0
|
||
for j in range(max_value):
|
||
if j <= i:
|
||
wA += hist[j]
|
||
uAtmp += j * hist[j]
|
||
else:
|
||
wB += hist[j]
|
||
uBtmp += j * hist[j]
|
||
if wA == 0:
|
||
wA = 1e-10
|
||
if wB == 0:
|
||
wB = 1e-10
|
||
uA = uAtmp / wA
|
||
uB = uBtmp / wB
|
||
u = uAtmp + uBtmp
|
||
|
||
# 计算类间方差
|
||
deltaTmp = wA * ((uA - u)**2) + wB * ((uB - u)**2)
|
||
# 找出最大类间方差以及阈值
|
||
if deltaTmp > deltaMax:
|
||
deltaMax = deltaTmp
|
||
threshold = i
|
||
|
||
# 二值化
|
||
det_img = img.copy()
|
||
det_img[img > threshold] = foreground
|
||
det_img[img <= threshold] = background
|
||
det_img[np.where(data_water_mask == 0)] = background
|
||
return det_img
|
||
|
||
|
||
def find_overexposure_area(img_path, threhold=4095):
|
||
# 第一步通过某个像素的光谱找到信号最强的波段
|
||
|
||
# 根据上步所得的波段号检测过曝区域
|
||
pass
|
||
|
||
|
||
@timeit
|
||
def find_severe_glint_area(img_path, water_mask, glint_wave=750, output_path=None):
|
||
if output_path is None:
|
||
output_path = append2filename(img_path, "_severe_glint_area")
|
||
|
||
glint_band_number = find_band_number(glint_wave, img_path)
|
||
|
||
dataset = gdal.Open(img_path)
|
||
num_bands = dataset.RasterCount
|
||
im_geotrans = dataset.GetGeoTransform()
|
||
im_proj = dataset.GetProjection()
|
||
|
||
tmp = dataset.GetRasterBand(glint_band_number + 1) # 波段计数从1开始
|
||
band_flare = tmp.ReadAsArray().astype(np.int16)
|
||
band_flare_stretch = (band_flare / band_flare.max() * 255).astype(np.int32)
|
||
|
||
dataset_water_mask = gdal.Open(water_mask)
|
||
data_water_mask = dataset_water_mask.GetRasterBand(1).ReadAsArray()
|
||
del dataset_water_mask
|
||
|
||
# 不加1报错：IndexError: index 73 is out of bounds for axis 0 with size 73
|
||
flare_binary = otsu(band_flare_stretch, band_flare_stretch.max() + 1, data_water_mask)
|
||
|
||
write_bands(img_path, output_path, flare_binary)
|
||
|
||
del dataset
|
||
|
||
return output_path
|
||
|
||
|
||
# Press the green button in the gutter to run the script.
|
||
if __name__ == '__main__':
|
||
img_path = r"D:\PycharmProjects\0water_rlx\test_data\ref_mosaic_1m_bsq"
|
||
|
||
parser = argparse.ArgumentParser(description="此程序通过大律法分割图像，提取耀斑最严重的区域，输出的栅格和输入的影像具有相同的行列数。")
|
||
|
||
parser.add_argument('-i1', '--input', type=str, required=True, help='输入影像文件的路径')
|
||
parser.add_argument('-i2', '--input_water_mask', type=str, required=True, help='输入水域掩膜文件的路径')
|
||
parser.add_argument('-gw', '--glint_wave', type=float, default=750.0, help='用于提取耀斑严重区域的波段.')
|
||
parser.add_argument('-o', '--output', type=str, help='输出文件的路径')
|
||
|
||
parser.add_argument('-v', '--verbose', action='store_true', help='启用详细模式')
|
||
|
||
args = parser.parse_args()
|
||
|
||
find_severe_glint_area(args.input, args.input_water_mask, args.glint_wave, args.output)
|