Commit d7fbce6b authored by luroth's avatar luroth
Browse files

Standalone Point Projection Tool

parent 32a3a238
# Photogrammetry Pipeline
# Copyright (C) 2019 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/>.
import sys
import csv
import os
import json
import pandas as pd
import numpy as np
import Metashape
app = Metashape.app
doc = app.document
def generate_image_masks_with_DEM():
path_polygons = app.getOpenFileName("Geojson with plot shapes", "*.geojson")
path_masks_output = app.getExistingDirectory("Output directory for image masks")
zero_gis = app.getBool("Inverted y: Should the origin be bottom left (for use in GIS)?")
chunk = doc.chunk
print("Processing chunk " + chunk.label)
points_list = []
with open(path_polygons) as json_file:
feature_collection = json.load(json_file)
features = feature_collection['features']
print(len(features), 'features in geojson found, extraction points for projection')
for feature in features:
properties = feature['properties']
geometry = feature['geometry']
label = properties['plot_label']
for i, corner in enumerate(geometry['coordinates'][0]):
points_list.append({'world_x': corner[0],
'world_y': corner[1],
'plot_label': label,
'point': i
})
print(len(points_list), 'points prepared for projection')
elevation = chunk.elevation
points = [Metashape.Vector((float(row['world_x']),
float(row['world_y']),
elevation.altitude(
Metashape.Vector((
float(row['world_x']),
float(row['world_y'])
)))
)) for row in points_list]
plot_labels = [row['plot_label'] for row in points_list]
point_numbers = [row['point'] for row in points_list]
rows = []
for i, point in enumerate(points):
for camera in chunk.cameras:
ret = camera.project(chunk.transform.matrix.inv().mulp(chunk.crs.unproject(point)))
if ret:
image_x, image_y = ret
if 0 <= image_y < camera.sensor.height and \
0 <= image_x < camera.sensor.width:
photo = camera.label
sensor = chunk.label
world_x, world_y, world_z = point
plot_label = plot_labels[i]
point_number = point_numbers[i]
row = {"image_x":image_x,
"image_y":image_y,
"image":photo,
"sensor":sensor,
"world_x":world_x,
"world_y":world_y,
"world_z": world_z,
"plot_label":plot_label,
"point_number":point_number,
"type":"canopy"}
rows.append(row)
keys = rows[0].keys()
with open(path_masks_output + "/projected_points.csv", 'w', newline='') as output_file:
dict_writer = csv.DictWriter(output_file, keys)
dict_writer.writeheader()
dict_writer.writerows(rows)
chunk.exportCameras(path_masks_output + "/camera_positions_.csv", format=Metashape.CamerasFormat.CamerasFormatOPK)
with open(path_masks_output + "/camera_positions_.csv", 'r') as infile:
with open(path_masks_output + "/camera_positions.csv", 'w') as outfile:
csvreader = csv.reader(infile, delimiter='\t', quotechar='|')
next(csvreader, None)
next(csvreader, None)
csvwriter = csv.writer(outfile, delimiter=',', quotechar='|', quoting=csv.QUOTE_MINIMAL)
csvwriter.writerow(
["image", "X", "Y", "Z", "Omega", "Phi", "Kappa", "r11", "r12", "r13", "r21", "r22", "r23", "r31",
"r32", "r33"])
csvwriter.writerows(row for row in csvreader)
# read projected points file
single_images_coords = pd.read_csv(path_masks_output + "/projected_points.csv", skip_blank_lines=True)
single_images_coords['image_y'] = -single_images_coords['image_y'] if zero_gis else single_images_coords['image_y']
# Group by polygons for local coords
single_images_coords_image = pd.DataFrame(single_images_coords.pivot_table(
index=['image', 'plot_label', 'type'],
values=['image_x', 'image_y'],
columns='point_number').to_records())
single_images_coords_image = single_images_coords_image.dropna()
# Group by images for world coords (mean)
single_images_coords_world = pd.DataFrame(single_images_coords.pivot_table(
index=['plot_label', 'type'],
values=['world_x', 'world_y', 'world_z'],
aggfunc="mean").to_records())
# Merge all to one df
single_images_coords_both = pd.merge(single_images_coords_image, single_images_coords_world,
on=["plot_label", "type"])
# Get camera position file to calculate viewpoints
camera_coords = pd.read_csv(path_masks_output + "/camera_positions.csv",)
# Merge camera position with projected plots
single_images_coords_all = pd.merge(single_images_coords_both, camera_coords, on="image")
# Calculate viewpoints
azimuth_proj = np.arctan2(single_images_coords_all.world_y - single_images_coords_all.Y,
single_images_coords_all.world_x - single_images_coords_all.X)
single_images_coords_all['azimuth_angle'] = 2 * np.pi - (azimuth_proj) % (2 * np.pi)
single_images_coords_all['zenith_angle'] = -(
np.arctan((single_images_coords_all.world_z - single_images_coords_all.Z) / np.sqrt(
(single_images_coords_all.world_x - single_images_coords_all.X) * (
single_images_coords_all.world_x - single_images_coords_all.X) +
(single_images_coords_all.world_y - single_images_coords_all.Y) *
(single_images_coords_all.world_y - single_images_coords_all.Y))) - np.pi / 2)
# Create polygons for geojson files
polygons = []
for index, row in single_images_coords_all.iterrows():
polygons.append(
[(row["('image_x', 0)"], row["('image_y', 0)"]),
(row["('image_x', 1)"], row["('image_y', 1)"]),
(row["('image_x', 2)"], row["('image_y', 2)"]),
(row["('image_x', 3)"], row["('image_y', 3)"]),
(row["('image_x', 0)"], row["('image_y', 0)"])]
)
# Geopandas df to allow geojson export
df_single_images_coords = pd.DataFrame(
single_images_coords_all[
['image', 'plot_label', 'type', 'world_x', 'world_y', 'world_z', 'azimuth_angle', 'zenith_angle']])
df_single_images_coords['geometry'] = polygons
# Group by image and export files
df_single_images_coords['image'] = df_single_images_coords['image'].str.split(".").apply(lambda x: x[0])
df_grouped = df_single_images_coords.groupby('image')
for image, df_group in df_grouped:
feature_collection = {}
feature_collection['type'] = "FeatureCollection"
features_list = []
dicts = df_group.to_dict('records')
for row in dicts:
properties = row
geometry = {"type": "Polygon", "coordinates": [row['geometry']]}
del properties['geometry']
feature = {"type": "Feature"}
feature['properties'] = properties
feature['geometry'] = geometry
features_list.append(feature)
feature_collection['features'] = features_list
with open(path_masks_output + "/" + image + ".geojson", 'w') as outfile:
json.dump(feature_collection, outfile, indent=4)
os.remove(path_masks_output + "/projected_points.csv")
os.remove(path_masks_output + "/camera_positions_.csv")
os.remove(path_masks_output + "/camera_positions.csv")
# Add menu item
app.addMenuItem("ETHZ CS Plugin/Generate image masks", generate_image_masks_with_DEM)
## Data Preparation
Basic requirements are
- Agisoft project with **aligned cameras** and a generated **DEM** in workspace
- Plot polygon file as **geojson** with the attribute `plot_label`
- Both Agisoft project and geojson need to be in the same *metric* coordinate system
Sample Agisoft project:
![Agisoft project](doc/agisoft_project.png)
Sample geojson:
```python
{"type": "FeatureCollection",
"features": [{
"type": "Feature",
"properties": {"plot_label": "FPWW0240067"},
"geometry": {
"type": "Polygon",
"coordinates": [
[[2693786.207215171, 1256288.427557322],
[2693790.4785625925, 1256285.030282317],
[2693791.407652783, 1256286.1984147034],
[2693787.1363053625, 1256289.5956897107],
[2693786.207215171, 1256288.427557322]]]
}
},
...
}
```
## Agisoft Metascan Preparation
You need to add pandas and numpy to the Agisoft Python installation. To do so, run the following line
in a windows comand window:
`"C:\Program Files\Agisoft\Metashape Pro\python\python.exe" -m pip install pandas==0.23 numpy`
## Generate image masks
1. Download 3. Select [Agisoft_standalone_project_points.py](Agisoft_standalone_project_points.py)
2. Open Agisoft Metascan
3. Select "Tools" -> "Run Script..." and run the downloaded script
4. The menu "ETHZ CS Plugin" appears.
5. Select "ETHZ CS Plugin" -> "Generate image masks"
6. Select the geojson file with the plot polygons that you prepared
7. Select an empty folder where the masks will be saved
8. Decide whether to use image coordinates (origin is in top left corner, usefull for processing in e.g. Python) or in bottom left corner (usefull for processing in GIS, e.g. QGIS)
If choosing reverse y coordinates, you may check the accuracy of the image mask in e.g. QGIS by drag-and-drop both origin image and generated plot mask:
![Maks](doc/mask_gis_style.png)
## Plot mask projection
[ImageProjectionAgisoft](ImageProjectionAgisoft/README.md)
\ No newline at end of file
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