radar-g/down.py

197 lines
7.9 KiB
Python

def _downsample_2d(src, mask, use_mask, method, fill_value, mode_rank, out):
src_w = src.shape[-1]
src_h = src.shape[-2]
out_w = out.shape[-1]
out_h = out.shape[-2]
if src_w == out_w and src_h == out_h:
return src
if out_w > src_w or out_h > src_h:
raise ValueError("invalid target size")
scale_x = src_w / out_w
scale_y = src_h / out_h
if method == DS_FIRST or method == DS_LAST:
for out_y in range(out_h):
src_yf0 = scale_y * out_y
src_yf1 = src_yf0 + scale_y
src_y0 = int(src_yf0)
src_y1 = int(src_yf1)
if src_y1 == src_yf1 and src_y1 > src_y0:
src_y1 -= 1
for out_x in range(out_w):
src_xf0 = scale_x * out_x
src_xf1 = src_xf0 + scale_x
src_x0 = int(src_xf0)
src_x1 = int(src_xf1)
if src_x1 == src_xf1 and src_x1 > src_x0:
src_x1 -= 1
done = False
value = fill_value
for src_y in range(src_y0, src_y1 + 1):
for src_x in range(src_x0, src_x1 + 1):
v = src[src_y, src_x]
if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
value = v
if method == DS_FIRST:
done = True
break
if done:
break
out[out_y, out_x] = value
elif method == DS_MODE:
max_value_count = int(scale_x + 1) * int(scale_y + 1)
values = np.zeros((max_value_count,), dtype=src.dtype)
frequencies = np.zeros((max_value_count,), dtype=np.uint32)
for out_y in range(out_h):
src_yf0 = scale_y * out_y
src_yf1 = src_yf0 + scale_y
src_y0 = int(src_yf0)
src_y1 = int(src_yf1)
wy0 = 1.0 - (src_yf0 - src_y0)
wy1 = src_yf1 - src_y1
if wy1 < _EPS:
wy1 = 1.0
if src_y1 > src_y0:
src_y1 -= 1
for out_x in range(out_w):
src_xf0 = scale_x * out_x
src_xf1 = src_xf0 + scale_x
src_x0 = int(src_xf0)
src_x1 = int(src_xf1)
wx0 = 1.0 - (src_xf0 - src_x0)
wx1 = src_xf1 - src_x1
if wx1 < _EPS:
wx1 = 1.0
if src_x1 > src_x0:
src_x1 -= 1
value_count = 0
for src_y in range(src_y0, src_y1 + 1):
wy = wy0 if (src_y == src_y0) else wy1 if (src_y == src_y1) else 1.0
for src_x in range(src_x0, src_x1 + 1):
wx = wx0 if (src_x == src_x0) else wx1 if (src_x == src_x1) else 1.0
v = src[src_y, src_x]
if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
w = wx * wy
found = False
for i in range(value_count):
if v == values[i]:
frequencies[i] += w
found = True
break
if not found:
values[value_count] = v
frequencies[value_count] = w
value_count += 1
w_max = -1.
value = fill_value
if mode_rank == 1:
for i in range(value_count):
w = frequencies[i]
if w > w_max:
w_max = w
value = values[i]
elif mode_rank <= max_value_count:
max_frequencies = np.full(mode_rank, -1.0, dtype=np.float64)
indices = np.zeros(mode_rank, dtype=np.int64)
for i in range(value_count):
w = frequencies[i]
for j in range(mode_rank):
if w > max_frequencies[j]:
max_frequencies[j] = w
indices[j] = i
break
value = values[indices[mode_rank - 1]]
out[out_y, out_x] = value
elif method == DS_MEAN:
for out_y in range(out_h):
src_yf0 = scale_y * out_y
src_yf1 = src_yf0 + scale_y
src_y0 = int(src_yf0)
src_y1 = int(src_yf1)
wy0 = 1.0 - (src_yf0 - src_y0)
wy1 = src_yf1 - src_y1
if wy1 < _EPS:
wy1 = 1.0
if src_y1 > src_y0:
src_y1 -= 1
for out_x in range(out_w):
src_xf0 = scale_x * out_x
src_xf1 = src_xf0 + scale_x
src_x0 = int(src_xf0)
src_x1 = int(src_xf1)
wx0 = 1.0 - (src_xf0 - src_x0)
wx1 = src_xf1 - src_x1
if wx1 < _EPS:
wx1 = 1.0
if src_x1 > src_x0:
src_x1 -= 1
v_sum = 0.0
w_sum = 0.0
for src_y in range(src_y0, src_y1 + 1):
wy = wy0 if (src_y == src_y0) else wy1 if (src_y == src_y1) else 1.0
for src_x in range(src_x0, src_x1 + 1):
wx = wx0 if (src_x == src_x0) else wx1 if (src_x == src_x1) else 1.0
v = src[src_y, src_x]
if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
w = wx * wy
v_sum += w * v
w_sum += w
if w_sum < _EPS:
out[out_y, out_x] = fill_value
else:
out[out_y, out_x] = v_sum / w_sum
elif method == DS_VAR or method == DS_STD:
for out_y in range(out_h):
src_yf0 = scale_y * out_y
src_yf1 = src_yf0 + scale_y
src_y0 = int(src_yf0)
src_y1 = int(src_yf1)
wy0 = 1.0 - (src_yf0 - src_y0)
wy1 = src_yf1 - src_y1
if wy1 < _EPS:
wy1 = 1.0
if src_y1 > src_y0:
src_y1 -= 1
for out_x in range(out_w):
src_xf0 = scale_x * out_x
src_xf1 = src_xf0 + scale_x
src_x0 = int(src_xf0)
src_x1 = int(src_xf1)
wx0 = 1.0 - (src_xf0 - src_x0)
wx1 = src_xf1 - src_x1
if wx1 < _EPS:
wx1 = 1.0
if src_x1 > src_x0:
src_x1 -= 1
v_sum = 0.0
w_sum = 0.0
wv_sum = 0.0
wvv_sum = 0.0
for src_y in range(src_y0, src_y1 + 1):
wy = wy0 if (src_y == src_y0) else wy1 if (src_y == src_y1) else 1.0
for src_x in range(src_x0, src_x1 + 1):
wx = wx0 if (src_x == src_x0) else wx1 if (src_x == src_x1) else 1.0
v = src[src_y, src_x]
if np.isfinite(v) and not (use_mask and mask[src_y, src_x]):
w = wx * wy
v_sum += v
w_sum += w
wv_sum += w * v
wvv_sum += w * v * v
if w_sum < _EPS:
out[out_y, out_x] = fill_value
else:
out[out_y, out_x] = (wvv_sum * w_sum - wv_sum * wv_sum) / w_sum / w_sum
if method == DS_STD:
out = np.sqrt(out).astype(out.dtype)
else:
raise ValueError('invalid upsampling method')
return out