Commit cdc797b4 authored by luroth's avatar luroth
Browse files

end of day

parent 045bfea0
...@@ -91,7 +91,7 @@ def point_in_poly(x,y,poly): ...@@ -91,7 +91,7 @@ def point_in_poly(x,y,poly):
return inside return inside
def ray_trace(line_vector, line_point, plane_point, plane_normal): def ray_trace(line_vector, line_point, plane_point, plane_normal):
epsilon=1e-6 epsilon=1e-8
ndotu = plane_normal.dot(line_vector) ndotu = plane_normal.dot(line_vector)
if abs(ndotu) < epsilon: if abs(ndotu) < epsilon:
...@@ -126,6 +126,9 @@ def intersect_line_with_planes(plane_normales, plane_points, line_vector, line_p ...@@ -126,6 +126,9 @@ def intersect_line_with_planes(plane_normales, plane_points, line_vector, line_p
# for all planes # for all planes
for plane_index in index_list: for plane_index in index_list:
if(plane_delta_dist[plane_index]>5):
print("?", end="\n" ,flush=True)
return [np.nan, np.nan, np.nan]
plane_normal = plane_normales[plane_index] plane_normal = plane_normales[plane_index]
plane_point = plane_points[plane_index, [0,1,2]] plane_point = plane_points[plane_index, [0,1,2]]
...@@ -136,17 +139,13 @@ def intersect_line_with_planes(plane_normales, plane_points, line_vector, line_p ...@@ -136,17 +139,13 @@ def intersect_line_with_planes(plane_normales, plane_points, line_vector, line_p
p1_x, p1_y, p2_x, p2_y, p3_x, p3_y = plane_points[plane_index, [0,1,3,4,6,7]].tolist() p1_x, p1_y, p2_x, p2_y, p3_x, p3_y = plane_points[plane_index, [0,1,3,4,6,7]].tolist()
hit = point_in_poly(intersect_point[0], intersect_point[1], [(p1_x, p1_y), (p2_x, p2_y), (p3_x, p3_y)] ) hit = point_in_poly(intersect_point[0], intersect_point[1], [(p1_x, p1_y), (p2_x, p2_y), (p3_x, p3_y)] )
if hit: if hit:
#print(".", end="\n") print('.', flush=True)
#print(".", end="", flush=True)
return intersect_point.tolist() return intersect_point.tolist()
else: else:
#print("_", end="") print("_", end="")
pass pass
print("?", end="\n" ,flush=True)
print("?")
return [np.nan, np.nan, np.nan] return [np.nan, np.nan, np.nan]
def main(argv): def main(argv):
...@@ -216,13 +215,13 @@ def main(argv): ...@@ -216,13 +215,13 @@ def main(argv):
# Coordinates of sample pixels of sensor (later used for interpolation), relative to center # Coordinates of sample pixels of sensor (later used for interpolation), relative to center
sensor_pixels_x = np.ceil(np.linspace(0, sensor_resolution_x-1, no_of_interpolation_points_x, dtype='int')) sensor_pixels_x = np.ceil(np.linspace(0, sensor_resolution_x-1, no_of_interpolation_points_x, dtype='int'))
sensor_pixels_y = np.ceil(np.linspace(0, sensor_resolution_y-1, no_of_interpolation_points_y, dtype='int')) sensor_pixels_y = np.ceil(np.linspace(0, sensor_resolution_y-1, no_of_interpolation_points_y, dtype='int'))
sensor_pixels_xy = np.transpose([np.tile(sensor_pixels_x, len(sensor_pixels_y)), np.repeat(sensor_pixels_y, len(sensor_pixels_x))]) sensor_pixels_xy = np.transpose([np.repeat(sensor_pixels_x, len(sensor_pixels_y)), np.tile(sensor_pixels_y, len(sensor_pixels_x))])
corner_index = [0, (no_of_interpolation_points_x-1), (no_of_interpolation_points_x*no_of_interpolation_points_y-1), (no_of_interpolation_points_x*no_of_interpolation_points_y-no_of_interpolation_points_x)] corner_index = [0, (no_of_interpolation_points_y-1), (no_of_interpolation_points_x*no_of_interpolation_points_y-1), (no_of_interpolation_points_x*no_of_interpolation_points_y-no_of_interpolation_points_y)]
sensor_pixels_coords_x = np.linspace((sensor_dimension_x / 2000), -(sensor_dimension_x / 2000), no_of_interpolation_points_x, dtype=np.float64) sensor_pixels_coords_x = (((sensor_pixels_x / sensor_resolution_x) - 0.5 ) * sensor_dimension_x) / 1000
sensor_pixels_coords_y = np.linspace((sensor_dimension_y / 2000), -(sensor_dimension_y / 2000), no_of_interpolation_points_y, dtype=np.float64) sensor_pixels_coords_y = (((sensor_pixels_y / sensor_resolution_y) - 0.5 ) * sensor_dimension_y) / 1000
sensor_pixels_coords_xy = np.transpose([np.tile(sensor_pixels_coords_x, len(sensor_pixels_coords_y)), np.repeat(sensor_pixels_coords_y, len(sensor_pixels_coords_x))]) sensor_pixels_coords_xy = np.transpose([np.repeat(sensor_pixels_coords_x, len(sensor_pixels_coords_y)), np.tile(sensor_pixels_coords_y, len(sensor_pixels_coords_x))])
sensor_pixel_coords_z = np.full((sensor_pixels_coords_xy.shape[0], 1), sensor_Z, dtype=np.float64) sensor_pixel_coords_z = np.full((sensor_pixels_coords_xy.shape[0], 1), sensor_Z, dtype=np.float64)
sensor_pixels_coords_xyz = np.concatenate((sensor_pixels_coords_xy,sensor_pixel_coords_z), axis=1) sensor_pixels_coords_xyz = np.concatenate((sensor_pixels_coords_xy,sensor_pixel_coords_z), axis=1)
...@@ -253,10 +252,6 @@ def main(argv): ...@@ -253,10 +252,6 @@ def main(argv):
focal_point = np.array([X, Y, Z]) focal_point = np.array([X, Y, Z])
print("rotate camera") print("rotate camera")
#rot_matrix = rotation_matrix(z_axis, kappa_rad).dot(rotation_matrix(y_axis, phi_rad)).dot(rotation_matrix(x_axis, omega_rad))
#rot_matrix = rotation_matrix(x_axis, omega_rad).dot(rotation_matrix(y_axis, phi_rad)).dot(rotation_matrix(z_axis, kappa_rad))
#sensor_pixels_vectors = np.dot(sensor_pixels_coords_xyz, rot_matrix)
sensor_pixels_vectors = np.array([np.dot(rotation_matrix(z_axis, kappa_rad), sensor_corner) for sensor_corner in sensor_pixels_coords_xyz]) sensor_pixels_vectors = np.array([np.dot(rotation_matrix(z_axis, kappa_rad), sensor_corner) for sensor_corner in sensor_pixels_coords_xyz])
sensor_pixels_vectors = np.array([np.dot(rotation_matrix(y_axis, phi_rad), sensor_corner) for sensor_corner in sensor_pixels_vectors]) sensor_pixels_vectors = np.array([np.dot(rotation_matrix(y_axis, phi_rad), sensor_corner) for sensor_corner in sensor_pixels_vectors])
sensor_pixels_vectors = np.array([np.dot(rotation_matrix(x_axis, omega_rad), sensor_corner) for sensor_corner in sensor_pixels_vectors]) sensor_pixels_vectors = np.array([np.dot(rotation_matrix(x_axis, omega_rad), sensor_corner) for sensor_corner in sensor_pixels_vectors])
...@@ -280,9 +275,9 @@ def main(argv): ...@@ -280,9 +275,9 @@ def main(argv):
shp.field('label') shp.field('label')
shp.poly(parts=[corner_polygon.tolist()], shapeType=shapefile.POLYGON) shp.poly(parts=[corner_polygon.tolist()], shapeType=shapefile.POLYGON)
shp.record('field_of_view') shp.record('field_of_view')
shp.save(os.path.join(output_directory, 'corners_' + PhotoID)) shp.save(os.path.join(output_directory, PhotoID + '_corners'))
prj = open(os.path.join(output_directory, 'corners_' + PhotoID + '.prj'), "w") prj = open(os.path.join(output_directory, PhotoID + '_corners.prj'), "w")
prj.write(epsg) prj.write(epsg)
prj.close() prj.close()
...@@ -304,57 +299,67 @@ def main(argv): ...@@ -304,57 +299,67 @@ def main(argv):
# interpolating # interpolating
print("interpolating") print("interpolating")
grid_y, grid_x = np.mgrid[0:sensor_resolution_y, 0:sensor_resolution_x] grid_x, grid_y = np.mgrid[0:sensor_resolution_x, 0:sensor_resolution_y]
sensor_pixels_xy_sel = sensor_pixels_xy[np.logical_not(np.isnan(ground_pixel_coords[:,0]))]
ground_pixel_coords_x = ground_pixel_coords[np.logical_not(np.isnan(ground_pixel_coords[:,0])),0]
ground_pixel_coords_x_interpol = griddata(sensor_pixels_xy_sel, ground_pixel_coords_x, (grid_x, grid_y), method='linear')
ground_pixel_coords_x = ground_pixel_coords[:,0] sensor_pixels_xy_sel = sensor_pixels_xy[np.logical_not(np.isnan(ground_pixel_coords[:,1]))]
ground_pixel_coords_x_interpol = griddata(sensor_pixels_xy, ground_pixel_coords_x, (grid_x, grid_y), method='linear') ground_pixel_coords_y = ground_pixel_coords[np.logical_not(np.isnan(ground_pixel_coords[:,1])),1]
ground_pixel_coords_y_interpol = griddata(sensor_pixels_xy_sel, ground_pixel_coords_y, (grid_x, grid_y), method='linear')
ground_pixel_coords_y = ground_pixel_coords[:,1] sensor_pixels_xy_sel = sensor_pixels_xy[np.logical_not(np.isnan(ground_pixel_coords[:,2]))]
ground_pixel_coords_y_interpol = griddata(sensor_pixels_xy, ground_pixel_coords_y, (grid_x, grid_y), method='linear') ground_pixel_coords_z = ground_pixel_coords[np.logical_not(np.isnan(ground_pixel_coords[:,2])),2]
ground_pixel_coords_z_interpol = griddata(sensor_pixels_xy_sel, ground_pixel_coords_z, (grid_x, grid_y), method='linear')
ground_pixel_coords_z = ground_pixel_coords[:,2] sensor_pixels_xy_sel = sensor_pixels_xy[np.logical_not(np.isnan(ground_pixel_angle_xy))]
ground_pixel_coords_z_interpol = griddata(sensor_pixels_xy, ground_pixel_coords_z, (grid_x, grid_y), method='linear') ground_pixel_angle_xy = ground_pixel_angle_xy[np.logical_not(np.isnan(ground_pixel_angle_xy))]
ground_pixel_angle_xy_interpol = griddata(sensor_pixels_xy_sel, ground_pixel_angle_xy, (grid_x, grid_y), method='linear')
ground_pixel_angle_xy_interpol = griddata(sensor_pixels_xy, ground_pixel_angle_xy, (grid_x, grid_y), method='linear') sensor_pixels_xy_sel = sensor_pixels_xy[np.logical_not(np.isnan(ground_pixel_angle_z))]
ground_pixel_angle_z_interpol = griddata(sensor_pixels_xy, ground_pixel_angle_z, (grid_x, grid_y), method='linear') ground_pixel_angle_z = ground_pixel_angle_z[np.logical_not(np.isnan(ground_pixel_angle_z))]
ground_pixel_angle_z_interpol = griddata(sensor_pixels_xy_sel, ground_pixel_angle_z, (grid_x, grid_y), method='linear')
print("create overview file") print("create overview file")
plt.figure(figsize=(10,7)) plt.figure(figsize=(10,7))
plt.subplot(231) plt.subplot(231)
plt.imshow(ground_pixel_coords_x_interpol) plt.imshow(ground_pixel_coords_x_interpol.T)
plt.title('X coord') plt.title('X coord')
plt.colorbar(orientation='horizontal') plt.colorbar(orientation='horizontal')
plt.subplot(232) plt.subplot(232)
plt.imshow(ground_pixel_coords_y_interpol) plt.imshow(ground_pixel_coords_y_interpol.T)
plt.title('Y coord') plt.title('Y coord')
plt.colorbar(orientation='horizontal') plt.colorbar(orientation='horizontal')
plt.subplot(233) plt.subplot(233)
plt.imshow(ground_pixel_coords_z_interpol) plt.imshow(ground_pixel_coords_z_interpol.T)
plt.title('Z coord') plt.title('Z coord')
plt.colorbar(orientation='horizontal') plt.colorbar(orientation='horizontal')
plt.subplot(234) plt.subplot(234)
plt.imshow(ground_pixel_angle_xy_interpol) plt.imshow(ground_pixel_angle_xy_interpol.T)
plt.title('xy view angle') plt.title('xy view angle')
plt.colorbar(orientation='horizontal') plt.colorbar(orientation='horizontal')
plt.subplot(235) plt.subplot(235)
plt.imshow(ground_pixel_angle_z_interpol) plt.imshow(ground_pixel_angle_z_interpol.T)
plt.title('z view angle') plt.title('z view angle')
plt.colorbar(orientation='horizontal') plt.colorbar(orientation='horizontal')
plt.savefig(os.path.join(output_directory, 'overview_' + PhotoID + '.png'), dpi = 80) plt.savefig(os.path.join(output_directory, PhotoID + '_overview.png'), dpi = 80)
plt.close()
print("\nwrite pixel matrix to TIFF") print("\nwrite pixel matrix to TIFF")
gdalnumeric.SaveArray(ground_pixel_coords_x_interpol, os.path.join(output_directory, 'metadata_coord_x_' + PhotoID + '.tiff'), format='GTiff' ) gdalnumeric.SaveArray(ground_pixel_coords_x_interpol.T, os.path.join(output_directory, PhotoID + '_metadata_coord_x.tif'), format='GTiff' )
gdalnumeric.SaveArray(ground_pixel_coords_y_interpol, os.path.join(output_directory, 'metadata_coord_y_' + PhotoID + '.tiff'), format='GTiff' ) gdalnumeric.SaveArray(ground_pixel_coords_y_interpol.T, os.path.join(output_directory, PhotoID + '_metadata_coord_y.tif'), format='GTiff' )
gdalnumeric.SaveArray(ground_pixel_coords_z_interpol, os.path.join(output_directory, 'metadata_coord_z_' + PhotoID + '.tiff'), format='GTiff' ) gdalnumeric.SaveArray(ground_pixel_coords_z_interpol.T, os.path.join(output_directory, PhotoID + '_metadata_coord_z.tif'), format='GTiff' )
gdalnumeric.SaveArray(ground_pixel_angle_xy_interpol, os.path.join(output_directory, 'metadata_angle_xy' + PhotoID + '.tiff'), format='GTiff' ) gdalnumeric.SaveArray(ground_pixel_angle_xy_interpol.T, os.path.join(output_directory, PhotoID + '_metadata_angle_xy.tif'), format='GTiff' )
gdalnumeric.SaveArray(ground_pixel_angle_z_interpol, os.path.join(output_directory, 'metadata_angle_z' + PhotoID + '.tiff'), format='GTiff' ) gdalnumeric.SaveArray(ground_pixel_angle_z_interpol.T, os.path.join(output_directory, PhotoID + '_metadata_angle_z.tif'), format='GTiff' )
print("\nDONE\n") print("\nDONE\n")
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment