To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 37e0fa25 authored by luroth's avatar luroth
Browse files

Merge branch 'master' into pre_release

# Conflicts:
#	Common/GISfunctions.py
#	README.md
parents 5dfb11e1 5dbd52c6
##############################################################################
#
# Common functions
# Project: PhenoFly UAS data processing tools
# File: Common functions
#
# Author: Lukas Roth (lukas.roth@usys.ethz.ch)
#
# Copyright (C) 2018 ETH Zürich, Lukas Roth (lukas.roth@usys.ethz.ch)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
# Dependency
import math
import numpy as np
import stl
from osgeo import gdal
......@@ -46,7 +61,7 @@ def ray_plane_intersect(ray_point, ray_vector, plane_point, plane_normal, epsilo
return psi
def image_rgb_cut(point, polygon):
def point_in_polygon(point, polygon):
"""Test if point is inside poylgon. Using
:param point: Tuple with coordinates of point to test
......@@ -126,10 +141,10 @@ def create_nband_GeoTiff(file_path, arrays, data_type, keys=None):
driver = gdal.GetDriverByName('GTiff')
file = file_path
ds = driver.Create(file, arrays[0].shape[1], arrays[0].shape[0], len(arrays), data_type )
ds = driver.Create(file, arrays[0].shape[1], arrays[0].shape[0], len(arrays), data_type)
for i, array in enumerate(arrays):
band_ = ds.GetRasterBand(i+1)
band_ = ds.GetRasterBand(i + 1)
band_.WriteArray(array)
if keys is not None:
for key, value in keys[i].items():
......
##############################################################################
#
# Direct georeference images based on camera position and Waypoint_flight
# Project: PhenoFly UAS data processing tools
# File: Direct georeference images based on camera position and Waypoint_flight
#
# Author: Lukas Roth (lukas.roth@usys.ethz.ch)
#
# Copyright (C) 2018 ETH Zürich, Lukas Roth (lukas.roth@usys.ethz.ch)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
##############################################################################
# Input:
# - Camera positions
# - Digital terrain model (DTM) in STL format
# - Camera/sensor parameter (sensor dimension, focal length, sensor resolution)
# - Number of interpolation points in x and y
# - EPSG coordinate system (only cartesian systems supported)
#
# Output:
# - Longitude/Latitude coordinates TIFF per image
# - Zenith/Azimuth angle TIFF per image
# - geoJSON showing shapes of all images
#
##############################################################################
# Dependencies
import csv
......@@ -30,10 +47,6 @@ import time
import math
import geojson
import multiprocessing
import matplotlib
import scipy
import osgeo
import shapely
# Specific imports
from multiprocessing import Process, Manager
......@@ -45,7 +58,6 @@ from geojson import Polygon, Feature, FeatureCollection
from Common.GISfunctions import ray_plane_intersect, point_in_polygon, find_bounding_box_mesh, get_triangle_normale, \
rotation_matrix, create_nband_GeoTiff
################################################################################
# Helper functions
################################################################################
......@@ -58,6 +70,7 @@ z_axis = [0, 0, 1]
# Max value in an uint16
uint16_max = 65536
def ray_tracing_worker(work_queue, result_queue, plane_normales, plane_points):
""" Worker to process ray tracing for pixel rays
......@@ -111,7 +124,9 @@ def process_ray_tracing_job(pixel_nr, plane_normales, plane_points, ray_vector,
dtm_approx_triangle_point = bounding_box[0]
# Estimate intersection of ray with approximation Waypoint_flight triangles
# Intersect ray with two triangles
approximation_point = np.array(ray_plane_intersect(ray_point, ray_vector, dtm_approx_triangle_point, dtm_approx_triangle_normale), dtype=np.float64)
approximation_point = np.array(
ray_plane_intersect(ray_point, ray_vector, dtm_approx_triangle_point, dtm_approx_triangle_normale),
dtype=np.float64)
# Log time for approximation
if verbose: t_approximation = round(time.time() - t_start, 2)
# If no approximation point found then point is outside bounding box, return NAN
......@@ -122,13 +137,16 @@ def process_ray_tracing_job(pixel_nr, plane_normales, plane_points, ray_vector,
# Calculate distance of planes to approximated point
# Simple euclidian distance
delta_plane_to_intersect = np.array(plane_points[:, [0, 1]], dtype=np.float64) - approximation_point[[0, 1]]
euclidian_distance_plane_to_intersect = (delta_plane_to_intersect[:, 0] ** 2) + (delta_plane_to_intersect[:, 1] ** 2)
euclidian_distance_plane_to_intersect = (delta_plane_to_intersect[:, 0] ** 2) + (
delta_plane_to_intersect[:, 1] ** 2)
# Sort planes by distance
# Create partitioned index, x closest planes in first partition
closest_planes_index = np.argpartition(euclidian_distance_plane_to_intersect,partition_size-1)[0:partition_size]
closest_planes_index = np.argpartition(euclidian_distance_plane_to_intersect, partition_size - 1)[0:partition_size]
# Recalculate intersection based on real Waypoint_flight
approximation_point_recalc = ray_plane_intersect(ray_point, ray_vector, plane_points[closest_planes_index[0]][0:3], plane_normales[closest_planes_index[0]])
approximation_point = np.array(approximation_point_recalc, dtype=np.float64) if (approximation_point_recalc is not None) else approximation_point
approximation_point_recalc = ray_plane_intersect(ray_point, ray_vector, plane_points[closest_planes_index[0]][0:3],
plane_normales[closest_planes_index[0]])
approximation_point = np.array(approximation_point_recalc, dtype=np.float64) if (
approximation_point_recalc is not None) else approximation_point
delta_plane_to_intersect = np.array(plane_points[:, [0, 1]], dtype=np.float64) - approximation_point[[0, 1]]
euclidian_distance_plane_to_intersect = (delta_plane_to_intersect[:, 0] ** 2) + (
......@@ -138,22 +156,30 @@ def process_ray_tracing_job(pixel_nr, plane_normales, plane_points, ray_vector,
if verbose: t_add = round(time.time() - t_start - t_approximation, 2)
# Intersect closest x planes with ray
i, intersect = process_ray_planes_intersect(closest_planes_index, plane_normales, plane_points, ray_vector, ray_point)
i, intersect = process_ray_planes_intersect(closest_planes_index, plane_normales, plane_points, ray_vector,
ray_point)
if intersect is None:
# No intersect with x closest planes. Use other planes too, up to the limit
if verbose:
print("WARNING: Pixel {0} needs full run. Think about increasing partitioning size (currently {1})".format(pixel_nr, partition_size), flush=True)
print("WARNING: Pixel {0} needs full run. Think about increasing partitioning size (currently {1})".format(
pixel_nr, partition_size), flush=True)
other_planes_index = np.argsort(euclidian_distance_plane_to_intersect)[partition_size:limit]
i, intersect = process_ray_planes_intersect(other_planes_index, plane_normales, plane_points, ray_vector, ray_point)
i, intersect = process_ray_planes_intersect(other_planes_index, plane_normales, plane_points, ray_vector,
ray_point)
# Log timespan until found
if verbose:
t_found = round(time.time()- t_start - t_approximation -t_add, 2)
t_found = round(time.time() - t_start - t_approximation - t_add, 2)
# Analyze run
if intersect is not None:
# Intersect found, print log
if verbose:
print("Pixel {0} calculated, {1} loops. Times: Init normal: {2}, Sort: {3}, Intersect: {4}".format(pixel_nr, i, t_approximation, t_add, t_found), flush=True)
print("Pixel {0} calculated, {1} loops. Times: Init normal: {2}, Sort: {3}, Intersect: {4}".format(pixel_nr,
i,
t_approximation,
t_add,
t_found),
flush=True)
elif pixel_nr % 50 == 0:
print("{} done".format(pixel_nr), flush=True)
# Return result
......@@ -161,8 +187,13 @@ def process_ray_tracing_job(pixel_nr, plane_normales, plane_points, ray_vector,
else:
# No intersect found, print log
if verbose:
t_found = round(time.time()- t_start - t_approximation -t_add, 2)
print("Pixel {0} not found, {1} loops. Times: Init normal: {2}, Sort: {3}, Intersect: {4}".format(pixel_nr, i, t_approximation, t_add, t_found), flush=True)
t_found = round(time.time() - t_start - t_approximation - t_add, 2)
print(
"Pixel {0} not found, {1} loops. Times: Init normal: {2}, Sort: {3}, Intersect: {4}".format(pixel_nr, i,
t_approximation,
t_add,
t_found),
flush=True)
elif pixel_nr % 50 == 0:
print("{} done".format(pixel_nr), flush=True)
# Return NAN array
......@@ -184,14 +215,14 @@ def process_ray_planes_intersect(index_list, plane_normales, plane_points, ray_v
for i, plane_index in enumerate(index_list):
# Get plane normale and point
plane_normal = plane_normales[plane_index]
plane_point = plane_points[plane_index, [0,1,2]]
plane_point = plane_points[plane_index, [0, 1, 2]]
# Intersect line with plane
intersect_point = ray_plane_intersect(ray_point, ray_vector, plane_point, plane_normal)
if(intersect_point is not None):
if (intersect_point is not None):
# Test if found intersection is in plane or outside
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_polygon((intersect_point[0], intersect_point[1]), [(p1_x, p1_y), (p2_x, p2_y), (p3_x, p3_y)] )
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_polygon((intersect_point[0], intersect_point[1]), [(p1_x, p1_y), (p2_x, p2_y), (p3_x, p3_y)])
if hit:
# If intersection is part of plane: True intersection, skip loop and return coordinates
return i, intersect_point.tolist()
......@@ -280,11 +311,11 @@ def process(camera_position_file, DTM_stl_file, no_of_interpolation_points_x, no
[np.repeat(sensor_pixels_x, len(sensor_pixels_y)), np.tile(sensor_pixels_y, len(sensor_pixels_x))])
sensor_pixels_x_dist = sensor_pixels_xy[:, 0] * (
1 + (sensor_pixels_xy[:, 0] ** 2 + sensor_pixels_xy[:, 1] ** 2) * k1) + cx * (
sensor_dimension_x / sensor_resolution_x)
1 + (sensor_pixels_xy[:, 0] ** 2 + sensor_pixels_xy[:, 1] ** 2) * k1) + cx * (
sensor_dimension_x / sensor_resolution_x)
sensor_pixels_y_dist = sensor_pixels_xy[:, 1] * (
1 + (sensor_pixels_xy[:, 0] ** 2 + sensor_pixels_xy[:, 1] ** 2) * k1) - cy * (
sensor_dimension_y / sensor_resolution_y)
1 + (sensor_pixels_xy[:, 0] ** 2 + sensor_pixels_xy[:, 1] ** 2) * k1) - cy * (
sensor_dimension_y / sensor_resolution_y)
sensor_pixels_pox_xy_dist = np.stack((sensor_pixels_x_dist, sensor_pixels_y_dist), axis=1)
......@@ -458,7 +489,7 @@ def process(camera_position_file, DTM_stl_file, no_of_interpolation_points_x, no
ground_pixel_angle_azimuth = 2 * np.pi - (azimuth_proj) % (2 * np.pi)
ground_pixel_angle_zenith = -(np.arctan((Z - sensor_pixel_coords[:, 2]) / np.sqrt(
(X - sensor_pixel_coords[:, 0]) * (X - sensor_pixel_coords[:, 0]) + (Y - sensor_pixel_coords[:, 1]) * (
Y - sensor_pixel_coords[:, 1]))) - np.pi / 2)
Y - sensor_pixel_coords[:, 1]))) - np.pi / 2)
# Interpolating x, y, z coordinates and viewing geometry for all sensor pixels
print("Interpolating pixels")
......@@ -567,6 +598,7 @@ def process(camera_position_file, DTM_stl_file, no_of_interpolation_points_x, no
print("\nDONE\n")
################################################################################
# Main
################################################################################
......@@ -615,8 +647,10 @@ def main(argv):
output_directory,
epsg, sensor_dimension_x, sensor_dimension_y, focal_length, sensor_resolution_x, sensor_resolution_y,
flight_height, verbose)
if __name__ == "__main__":
main(sys.argv[1:])
print("\n================================================================================\nAll photos processed, END\n")
print(
"\n================================================================================\nAll photos processed, END\n")
##############################################################################
#
# Extract samples from direct georeferenced images using a samples mask
# Project: PhenoFly UAS data processing tools
# File: Extract samples from direct georeferenced images using a samples mask
#
# Author: Lukas Roth (lukas.roth@usys.ethz.ch)
#
# Copyright (C) 2018 ETH Zürich, Lukas Roth (lukas.roth@usys.ethz.ch)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
##############################################################################
#
# Input:
# - RAW images
......@@ -14,7 +33,6 @@
# - Sampled image files as georeferenced geoTIFF including 1bit mask
# - Zenith/Azimuth angle TIFF as georeferenced geoTIFF including 1bit mask
#
##############################################################################
# Dependencies
import sys
......@@ -30,9 +48,6 @@ from matplotlib import path
import geojson
import math
import pandas as pd
from shapely.geometry import Polygon
import gdalnumeric
from Common.GISfunctions import create_nband_GeoTiff
import osr
from affine import Affine
......
##############################################################################
#
# Coverts an Agisoft camera position file to a regular CSV file
# Project: PhenoFly UAS data processing tools
# File: Coverts an Agisoft camera position file to a regular CSV file
#
# Author: Lukas Roth (lukas.roth@usys.ethz.ch)
#
# Copyright (C) 2018 ETH Zürich, Lukas Roth (lukas.roth@usys.ethz.ch)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Based on:
# Roth, Aasen, Walter, Liebisch 2018: Extracting leaf area index using viewing
# geometry effects—A new perspective on high-resolution unmanned aerial system
# photography, ISPRS Journal of Photogrammetry and Remote Sensing.
# https://doi.org/10.1016/j.isprsjprs.2018.04.012
#
##############################################################################
# Input:
# - Agisoft camera position file
#
# Output:
# - XYZ Omega Phi Kappa CSV file
#
##############################################################################
import csv
......
##############################################################################
#
# Sample script to project iamges on DTM followed by plot-based sample extraction
# Project: PhenoFly UAS data processing tools
# File: Sample script to project iamges on DTM followed by plot-based sample extraction
#
# Author: Lukas Roth (lukas.roth@usys.ethz.ch)
#
# Copyright (C) 2018 ETH Zürich, Lukas Roth (lukas.roth@usys.ethz.ch)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Based on:
# Roth, Aasen, Walter, Liebisch 2018: Extracting leaf area index using viewing
# geometry effects—A new perspective on high-resolution unmanned aerial system
# photography, ISPRS Journal of Photogrammetry and Remote Sensing.
# https://doi.org/10.1016/j.isprsjprs.2018.04.012
#
##############################################################################
# Input:
# - Agisoft camera position file (_TestData/sonyDX100II/camera_position_RGB.txt)
# - Digital terrain model in STL file format (_TestData/sonyDX100II/DTM.stl)
......@@ -22,8 +40,6 @@
# - Extracted plot-based image parts (_Demo_output/samples/)
# - Viewpoint of each extracted sample (_Demo_output/samples/viewpoint.csv)
#
##############################################################################
from ImageProjectionSampling import ReformatAgisoftCameraPositionFile
from ImageProjectionSampling import DirectGeoreference
......
Markdown is supported
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