Skip to content
Snippets Groups Projects
Commit f791cc31 authored by Narayanarao Bhogapurapu's avatar Narayanarao Bhogapurapu
Browse files

mask update

parent 9f795912
No related branches found
No related tags found
No related merge requests found
File deleted
File deleted
File deleted
File deleted
File deleted
......@@ -248,6 +248,7 @@ def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_):
cohBlock = cohArray[start_i:end_i, start_j:end_j]
if np.count_nonzero(~np.isnan(cohBlock)) <5:
# print('no coh')
count = np.count_nonzero(~np.isnan(lidarBlock))
s_parm = 0
c_parm = 0
......@@ -259,22 +260,23 @@ def process_block(i, j, cohArray, lidarArray, initial_ws, htl, htg, parm_):
return start_i, end_i, start_j, end_j, s_parm, c_parm, rmse_parm, ht, count
# elif (np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 1) or max(lidarBlock.shape)>512:
elif (np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 1) or max(lidarBlock.shape)>512:
elif (np.isfinite(lidarBlock).any() and np.count_nonzero(~np.isnan(lidarBlock)) > 3) or max(lidarBlock.shape)>512:
parm = cal_(cohBlock, lidarBlock, htl, htg)
mask = np.where(~np.isnan(lidarBlock), 1, 0)
if np.all(np.array(parm) == 0):
parm = parm_.copy()
del parm_
# if np.all(np.array(parm) == 0):
# parm = parm_.copy()
# del parm_
if lidarBlock.shape[0]*lidarBlock.shape[1]>initial_ws*initial_ws:
# mask[np.shape(mask)[0]//2,np.shape(mask)[1]//2]=1
print('blank_blocks',lidarBlock.shape,parm)
np.fill_diagonal(mask, 1)
mask = np.flipud(mask)
np.fill_diagonal(mask, 1)
s_parm = np.full(lidarBlock.shape, parm[1]) * mask
c_parm = np.full(lidarBlock.shape, parm[2]) * mask
s_parm = np.full(lidarBlock.shape, parm[1])
c_parm = np.full(lidarBlock.shape, parm[2])
rmse_parm = np.full(lidarBlock.shape, parm[4])
count = np.count_nonzero(~np.isnan(lidarBlock))
......
......@@ -28,7 +28,7 @@ from sklearn.model_selection import train_test_split
from plotting import density_scatter, plot_data, rmse
from algo import arc_sinc_fast,arc_sinc, arc_sinc_, cal_, dynamicWindow
from rst_io import read_bin, write_tif, write_bin
from utils import blockshaped, unblockshaped, unpad, spatial_intp_lin, split_sample
from utils import blockshaped, unblockshaped, unpad, spatial_intp_lin, split_sample, spatial_intp_idw
# Initialize EMPatches instance
emp = EMPatches()
......@@ -148,6 +148,26 @@ def rvog_inverse(args):
t1 = time.time()
s, c, rmse__, ht_, count = dynamicWindow(cor, lidar_ht_cal, args.window_size, args.htl, args.htg)
temp_lidar = blockshaped(lidar_ht_cal, args.window_size, args.window_size)
temp_mask = np.zeros(temp_lidar.shape)
for win in tqdm(range(np.shape(temp_lidar)[0])):
mask = temp_lidar[win,:,:].copy()
mask[~np.isnan(mask)] = 1
temp_mask[win,:,:] = mask
if np.all(temp_lidar[win,:,:]==0) or np.all(np.isnan(temp_lidar[win,:,:])):
mask = np.zeros(temp_mask[win,:,:].shape)
np.fill_diagonal(mask, 1)
mask = np.flipud(mask)
np.fill_diagonal(mask, 1)
temp_mask[win,:,:] = mask
temp_mask = unblockshaped(temp_mask, rows,cols)
s = s*temp_mask
c = c*temp_mask
# Mask out invalid values
ht_[ht_ == 0] = np.nan
s[s == 0] = np.nan
......@@ -171,10 +191,18 @@ def rvog_inverse(args):
t2 = time.time()
print("[2/3] Generating height map... ",end=" ")
# Interpolate calibration parameters
# # Interpolate calibration parameters
ci = spatial_intp_lin(c)
si = spatial_intp_lin(s)
# # Interpolate calibration parameters
# ci = spatial_intp_idw(c)
# si = spatial_intp_idw(s)
# # Interpolate calibration parameters
# ci = spatial_intp_kriging(c)
# si = spatial_intp_kriging(s)
# Clip extreme values
si[si > 1.6] = 1.6
ci[ci > 16] = 16
......
......@@ -21,7 +21,7 @@ from scipy.interpolate import griddata
from tqdm import tqdm
from scipy.interpolate import interpn
from sklearn.model_selection import train_test_split
# from pykrige.ok import OrdinaryKriging
def blockshaped(arr, nrows, ncols):
"""
......@@ -73,6 +73,52 @@ def spatial_intp_lin(array):
return intp
# def spatial_intp_kriging(array):
# x = np.arange(0, array.shape[1])
# y = np.arange(0, array.shape[0])
# # Mask invalid values
# array = np.ma.masked_invalid(array)
# xx, yy = np.meshgrid(x, y)
# # Get only the valid values
# x1 = xx[~array.mask]
# y1 = yy[~array.mask]
# newarr = array[~array.mask]
# # Perform ordinary kriging
# OK = OrdinaryKriging(x1, y1, newarr, variogram_model='linear', verbose=False, enable_plotting=False)
# z, ss = OK.execute('grid', x, y)
# return z
def spatial_intp_idw(array, power=0.5):
x = np.arange(0, array.shape[1])
y = np.arange(0, array.shape[0])
# Mask invalid values
array = np.ma.masked_invalid(array)
xx, yy = np.meshgrid(x, y)
# Get only the valid values
x1 = xx[~array.mask]
y1 = yy[~array.mask]
newarr = array[~array.mask]
# Function to compute IDW
def idw(xi, yi, x1, y1, values, power):
dist = np.sqrt((x1 - xi) ** 2 + (y1 - yi) ** 2)
# Handle case where distance is zero
if np.any(dist == 0):
return values[dist == 0][0]
weights = 1 / (dist ** power)
return np.sum(weights * values) / np.sum(weights)
# Apply IDW for each point in the grid
idw_values = np.array([[idw(i, j, x1, y1, newarr, power) for i in x] for j in y])
return idw_values
def split_sample(a,frac=0.3):
"""
Returns calibration and validation 2d arrays
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment