Skip to content

Commit

Permalink
Support tsai camera loading and updates to gallery plotter (#74)
Browse files Browse the repository at this point in the history 
This commit adds the ability to properly load tsai camera files. There are also minor updates to fix some style errors (linting) across the repo, and to make the Gallery plotting script functional again - a bandaid until we update and improve it to include in the package.

---------

Co-authored-by: Ben Purinton <[email protected]>
  • Loading branch information
@dshean @bpurinton
dshean and bpurinton authored Feb 5, 2025
1 parent db6bd46 commit 1600d70
Show file tree
Hide file tree
Showing 2 changed files with 69 additions and 31 deletions.
75 changes: 52 additions & 23 deletions asp_plot/csm_camera.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -215,7 +215,7 @@ def format_stat_value(value):
return f"{value:.2e}" if abs(value) < 0.01 else f"{value:.2f}"


def plot_stats_text(ax, mean, std, unit='m'):
def plot_stats_text(ax, mean, std, unit="m"):
"""
Plots a text annotation on the given axis displaying the mean and standard deviation of a value.
Expand Down Expand Up @@ -449,7 +449,7 @@ def csm_camera_summary_plot(
lw=1,
label="X position (easting)",
)
plot_stats_text(ax1, cam1_x_position_diff_mean, cam1_x_position_diff_std, unit='m')
plot_stats_text(ax1, cam1_x_position_diff_mean, cam1_x_position_diff_std, unit="m")
ax2 = axes[0, 2]
ax2.plot(
frame_cam1,
Expand All @@ -458,7 +458,7 @@ def csm_camera_summary_plot(
lw=1,
label="Y position (northing)",
)
plot_stats_text(ax2, cam1_y_position_diff_mean, cam1_y_position_diff_std, unit='m')
plot_stats_text(ax2, cam1_y_position_diff_mean, cam1_y_position_diff_std, unit="m")
ax3 = axes[0, 3]
ax3.plot(
frame_cam1,
Expand All @@ -467,7 +467,7 @@ def csm_camera_summary_plot(
lw=1,
label="Z position (altitude)",
)
plot_stats_text(ax3, cam1_z_position_diff_mean, cam1_z_position_diff_std, unit='m')
plot_stats_text(ax3, cam1_z_position_diff_mean, cam1_z_position_diff_std, unit="m")

# Share y-axis for position diff plots
min_val_position_diff = min(
Expand Down Expand Up @@ -525,7 +525,7 @@ def csm_camera_summary_plot(
linestyle="--",
label="Original Roll",
)
plot_stats_text(ax1_r, cam1_roll_diff_mean, cam1_roll_diff_std, unit='°')
plot_stats_text(ax1_r, cam1_roll_diff_mean, cam1_roll_diff_std, unit="°")

ax2 = axes[1, 2]
ax2.plot(frame_cam1, gdf_cam1.pitch_diff, c="#FFA500", lw=1, label="Pitch Diff")
Expand All @@ -538,7 +538,7 @@ def csm_camera_summary_plot(
linestyle="--",
label="Original Pitch",
)
plot_stats_text(ax2_r, cam1_pitch_diff_mean, cam1_pitch_diff_std, unit='°')
plot_stats_text(ax2_r, cam1_pitch_diff_mean, cam1_pitch_diff_std, unit="°")

ax3 = axes[1, 3]
ax3.plot(frame_cam1, gdf_cam1.yaw_diff, c="#FFB347", lw=1, label="Yaw Diff")
Expand All @@ -551,7 +551,7 @@ def csm_camera_summary_plot(
linestyle="--",
label="Original Yaw",
)
plot_stats_text(ax3_r, cam1_yaw_diff_mean, cam1_yaw_diff_std, unit='°')
plot_stats_text(ax3_r, cam1_yaw_diff_mean, cam1_yaw_diff_std, unit="°")

# Share y-axis for angular diff plots
min_val_angle_diff = min(
Expand Down Expand Up @@ -661,7 +661,9 @@ def csm_camera_summary_plot(
lw=1,
label="X position (easting)",
)
plot_stats_text(ax1, cam2_x_position_diff_mean, cam2_x_position_diff_std, unit='m')
plot_stats_text(
ax1, cam2_x_position_diff_mean, cam2_x_position_diff_std, unit="m"
)
ax2 = axes[2, 2]
ax2.plot(
frame_cam2,
Expand All @@ -670,7 +672,9 @@ def csm_camera_summary_plot(
lw=1,
label="Y position (northing)",
)
plot_stats_text(ax2, cam2_y_position_diff_mean, cam2_y_position_diff_std, unit='m')
plot_stats_text(
ax2, cam2_y_position_diff_mean, cam2_y_position_diff_std, unit="m"
)
ax3 = axes[2, 3]
ax3.plot(
frame_cam2,
Expand All @@ -679,7 +683,9 @@ def csm_camera_summary_plot(
lw=1,
label="Z position (altitude)",
)
plot_stats_text(ax3, cam2_z_position_diff_mean, cam2_z_position_diff_std, unit='m')
plot_stats_text(
ax3, cam2_z_position_diff_mean, cam2_z_position_diff_std, unit="m"
)

# Share y-axis for position diff plots
min_val_position_diff = min(
Expand Down Expand Up @@ -739,7 +745,7 @@ def csm_camera_summary_plot(
linestyle="--",
label="Original Roll",
)
plot_stats_text(ax1_r, cam2_roll_diff_mean, cam2_roll_diff_std, unit='°')
plot_stats_text(ax1_r, cam2_roll_diff_mean, cam2_roll_diff_std, unit="°")

ax2 = axes[3, 2]
ax2.plot(frame_cam2, gdf_cam2.pitch_diff, c="#FFA500", lw=1, label="Pitch Diff")
Expand All @@ -752,7 +758,7 @@ def csm_camera_summary_plot(
linestyle="--",
label="Original Pitch",
)
plot_stats_text(ax2_r, cam2_pitch_diff_mean, cam2_pitch_diff_std, unit='°')
plot_stats_text(ax2_r, cam2_pitch_diff_mean, cam2_pitch_diff_std, unit="°")

ax3 = axes[3, 3]
ax3.plot(frame_cam2, gdf_cam2.yaw_diff, c="#FFB347", lw=1, label="Yaw Diff")
Expand All @@ -765,7 +771,7 @@ def csm_camera_summary_plot(
linestyle="--",
label="Original Yaw",
)
plot_stats_text(ax3_r, cam2_yaw_diff_mean, cam2_yaw_diff_std, unit='°')
plot_stats_text(ax3_r, cam2_yaw_diff_mean, cam2_yaw_diff_std, unit="°")

# Share y-axis for angular diff plots
min_val_angle_diff = min(
Expand Down Expand Up @@ -986,8 +992,8 @@ def read_csm_cam(json_file):

def read_tsai_cam(tsai):
"""
read tsai frame model from asp and return a python dictionary containing the parameters
See ASP's frame camera implementation here: https://stereopipeline.readthedocs.io/en/latest/pinholemodels.html
read tsai frame camera model and return a python dictionary containing the parameters
See ASP documentation: https://stereopipeline.readthedocs.io/en/latest/pinholemodels.html
Parameters
----------
tsai: str
Expand All @@ -1002,32 +1008,55 @@ def read_tsai_cam(tsai):
with open(tsai, "r") as f:
content = f.readlines()
content = [x.strip() for x in content]
fu = np.float64(content[2].split(" = ", 4)[1]) # focal length in x
fv = np.float64(content[3].split(" = ", 4)[1]) # focal length in y
cu = np.float64(content[4].split(" = ", 4)[1]) # optical center in x
cv = np.float64(content[5].split(" = ", 4)[1]) # optical center in y
fu = float(content[2].split(" = ", 4)[1]) # focal length in x
fv = float(content[3].split(" = ", 4)[1]) # focal length in y
cu = float(content[4].split(" = ", 4)[1]) # optical center in x
cv = float(content[5].split(" = ", 4)[1]) # optical center in y
cam = content[9].split(" = ", 10)[1].split(" ")
cam_cen = [np.float64(x) for x in cam] # camera center coordinates in ECEF
cam_cen = [float(x) for x in cam] # camera center coordinates in ECEF
rot = content[10].split(" = ", 10)[1].split(" ")
rot_mat = [
np.float64(x) for x in rot
float(x) for x in rot
] # rotation matrix for camera to world coordinates transformation

# Reshape as 3x3 matrix
rot_mat = np.reshape(rot_mat, (3, 3))

pitch = np.float64(content[11].split(" = ", 10)[1]) # pixel pitch
pitch = float(content[11].split(" = ", 10)[1]) # pixel pitch
tsai_dict = {
"camera": camera,
"focal_length": (fu, fv),
"optical_center": (cu, cv),
"cam_cen_ecef": cam_cen,
"cam_cen_ecef": Point(cam_cen),
"rotation_matrix": rot_mat,
"pitch": pitch,
}
return tsai_dict


def tsai_list_to_gdf(tsai_fn_list):
"""
Read a list of tsai camera files and return a GeoDataFrame
Parameters
----------
tsai_fn_list: list
List of tsai filenames
Returns
----------
output: GeoDataFrame
GeoPandas GeoDataFrame containing camera model parameters
"""
tsai_dict_list = []
for tsai_fn in tsai_fn_list:
tsai_dict = read_tsai_cam(tsai_fn)
tsai_dict_list.append(tsai_dict)

gdf = gpd.GeoDataFrame(tsai_dict_list, geometry="cam_cen_ecef", crs="EPSG:4978")
gdf.set_index("camera")

return gdf


def read_frame_csm_cam(json_file):
"""
Read rotation from a CSM Frame json state file.
Expand Down
25 changes: 17 additions & 8 deletions original_code/gallery.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
from pygeotools.lib import iolib, timelib, geolib, malib
from imview.lib import pltlib

add_cbar = False
add_cbar = True
hs = False

#For DEMs
Expand Down Expand Up @@ -50,10 +50,13 @@
#nrows = int(np.sqrt(n))+1
#ncols = nrows

#For projected images, set to True
shared_extent=False

if add_cbar:
grid = ImageGrid(f, 111, nrows_ncols=(nrows, ncols), axes_pad=(0.05,0.2), share_all=True, cbar_mode='single', cbar_pad=0.1, cbar_set_cax=False)
grid = ImageGrid(f, 111, nrows_ncols=(nrows, ncols), axes_pad=(0.05,0.2), share_all=shared_extent, cbar_mode='single', cbar_pad=0.1, cbar_set_cax=False)
else:
grid = ImageGrid(f, 111, nrows_ncols=(nrows, ncols), axes_pad=(0.05,0.2), share_all=True, cbar_mode=None)
grid = ImageGrid(f, 111, nrows_ncols=(nrows, ncols), axes_pad=(0.05,0.2), share_all=shared_extent, cbar_mode=None)

#f, axa = plt.subplots(nrows, ncols, sharex='all', sharey='all', figsize=(10,10))
#plt.subplots_adjust(wspace=0, hspace=0.2)
Expand All @@ -76,7 +79,7 @@
#dem_clim = (2934, 3983)
hs_clim = (1, 255)

for i,dem_fn in enumerate(dem_fn_list):
for i, dem_fn in enumerate(dem_fn_list):
ax = grid[i]
print(dem_fn)
dem_ds = iolib.fn_getds(dem_fn)
Expand All @@ -88,11 +91,17 @@
else:
dem_hs = geolib.gdaldem_mem_ds(dem_ds, 'hillshade', returnma=True)
hs_im = ax.imshow(dem_hs, vmin=hs_clim[0], vmax=hs_clim[1], cmap='gray')
dt = timelib.fn_getdatetime(dem_fn)
#dt = timelib.fn_getdatetime(dem_fn)
dt = None
if dt is not None:
title = dt.strftime('%Y-%m-%d')
t = ax.set_title(title, fontdict={'fontsize':6})
t.set_position([0.5, 0.95])
else:
#Manually adjust based on structure of filename
#title = os.path.splitext(dem_fn)[0]
title = os.path.splitext(dem_fn)[0].split('_')[0]
#title = ''.join(dem_fn.split('-')[:3])
t = ax.set_title(title, fontdict={'fontsize':5})
t.set_position([0.5, 0.95])
dem_im = ax.imshow(dem, vmin=dem_clim[0], vmax=dem_clim[1], cmap=cmap, alpha=alpha)
#ax.set_facecolor('k')
ax.set_facecolor('0.5')
Expand All @@ -109,7 +118,7 @@
cbar_kwargs = {'extend':'both', 'alpha':1.0}
cbar = grid.cbar_axes[0].colorbar(dem_im, **cbar_kwargs)
cbar.update_normal(dem_im)
cbar.set_label(cbar_lbl)
#cbar.set_label(cbar_lbl)

#res = geolib.get_res(dem_ds)[0]
#pltlib.add_scalebar(grid[-1], res=res)
Expand Down

0 comments on commit 1600d70

Please sign in to comment.