diff --git a/app/api/geoutils.py b/app/api/geoutils.py deleted file mode 100644 index ef337a0f..00000000 --- a/app/api/geoutils.py +++ /dev/null @@ -1,59 +0,0 @@ -import rasterio.warp -import numpy as np -from rasterio.crs import CRS - -# GEOS has some weird bug where -# we can't simply call geom.tranform(srid) -# so we write our own - -def geom_transform_wkt_bbox(geom, dataset, bbox_crs="geographic", wkt_crs="raster"): - """ - :param geom GEOSGeometry - :param dataset rasterio dataset - :param bbox_crs CRS of bbox (geographic --> lat/lon | projected --> north/east | raster --> pixels) - :param wkt_crs CRS of WKT (raster --> pixels | projected --> north/east) - :return (WKT, bbox) - """ - if not geom.srid: - raise ValueError("Geometry must have an SRID") - if not dataset.crs: - raise ValueError("Dataset must have a CRS") - - coords = geom.tuple - if len(coords) == 1: - xs, ys = zip(*coords[0]) - tx, ty = rasterio.warp.transform(CRS.from_epsg(geom.srid), dataset.crs, xs, ys) - raster_coords = [dataset.index(x, y, op=np.round) for x, y in zip(tx, ty)] - - if bbox_crs == 'geographic': - minx = min(xs) - maxx = max(xs) - miny = min(ys) - maxy = max(ys) - elif bbox_crs == 'projected': - minx = min(tx) - maxx = max(tx) - miny = min(ty) - maxy = max(ty) - elif bbox_crs == 'raster': - coords = np.array(raster_coords) - px = coords[:, 1] - py = coords[:, 0] - minx = px.min() - maxx = px.max() - miny = py.min() - maxy = py.max() - else: - raise ValueError("Invalid bbox_crs") - - if wkt_crs == "raster": - out = ", ".join(f"{x} {y}" for y, x in raster_coords) - elif wkt_crs == "projected": - out = ", ".join(f"{x} {y}" for x, y in zip(tx, ty)) - else: - raise ValueError("Invalid wkt_crs") - - wkt = f"POLYGON (({out}))" - return wkt, (minx, miny, maxx, maxy) - else: - raise ValueError("Cannot transform complex geometries to WKT") diff --git a/app/api/tasks.py b/app/api/tasks.py index 21380aa7..acd65030 100644 --- a/app/api/tasks.py +++ b/app/api/tasks.py @@ -36,7 +36,7 @@ from .tags import TagsField from app.security import path_traversal_check from django.utils.translation import gettext_lazy as _ from .fields import PolygonGeometryField -from .geoutils import geom_transform_wkt_bbox +from app.geoutils import geom_transform_wkt_bbox from webodm import settings def flatten_files(request_files): diff --git a/app/api/tiler.py b/app/api/tiler.py index 471e29ae..0776100d 100644 --- a/app/api/tiler.py +++ b/app/api/tiler.py @@ -25,7 +25,7 @@ from .hsvblend import hsv_blend from .hillshade import LightSource from .formulas import lookup_formula, get_algorithm_list, get_auto_bands from .tasks import TaskNestedView -from .geoutils import geom_transform_wkt_bbox +from app.geoutils import geom_transform_wkt_bbox from rest_framework import exceptions from rest_framework.response import Response from worker.tasks import export_raster, export_pointcloud diff --git a/app/geoutils.py b/app/geoutils.py new file mode 100644 index 00000000..816a8f2b --- /dev/null +++ b/app/geoutils.py @@ -0,0 +1,81 @@ +import rasterio.warp +import numpy as np +from rasterio.crs import CRS + +# GEOS has some weird bug where +# we can't simply call geom.tranform(srid) +# so we write our own + +def geom_transform_wkt_bbox(geom, dataset, bbox_crs="geographic", wkt_crs="raster"): + """ + :param geom GEOSGeometry + :param dataset rasterio dataset | path to raster + :param bbox_crs CRS of bbox (geographic --> lat/lon | projected --> north/east | raster --> pixels) + :param wkt_crs CRS of WKT (raster --> pixels | projected --> north/east) + :return (WKT, bbox) + """ + if not geom.srid: + raise ValueError("Geometry must have an SRID") + + close_ds = False + if isinstance(dataset, str): + dataset = rasterio.open(dataset, 'r') + close_ds = True + + try: + if not dataset.crs: + raise ValueError("Dataset must have a CRS") + + coords = geom.tuple + if len(coords) == 1: + xs, ys = zip(*coords[0]) + tx, ty = rasterio.warp.transform(CRS.from_epsg(geom.srid), dataset.crs, xs, ys) + raster_coords = [dataset.index(x, y, op=np.round) for x, y in zip(tx, ty)] + + if bbox_crs == 'geographic': + minx = min(xs) + maxx = max(xs) + miny = min(ys) + maxy = max(ys) + elif bbox_crs == 'projected': + minx = min(tx) + maxx = max(tx) + miny = min(ty) + maxy = max(ty) + elif bbox_crs == 'raster': + coords = np.array(raster_coords) + px = coords[:, 1] + py = coords[:, 0] + minx = px.min() + maxx = px.max() + miny = py.min() + maxy = py.max() + else: + raise ValueError("Invalid bbox_crs") + + if wkt_crs == "raster": + out = ", ".join(f"{x} {y}" for y, x in raster_coords) + elif wkt_crs == "projected": + out = ", ".join(f"{x} {y}" for x, y in zip(tx, ty)) + else: + raise ValueError("Invalid wkt_crs") + + wkt = f"POLYGON (({out}))" + return wkt, (minx, miny, maxx, maxy) + else: + raise ValueError("Cannot transform complex geometries to WKT") + finally: + if close_ds: + dataset.close() + +def geom_transform(geom, epsg): + if not geom.srid: + raise ValueError("Geometry must have an SRID") + + coords = geom.tuple + if len(coords) == 1: + xs, ys = zip(*coords[0]) + tx, ty = rasterio.warp.transform(CRS.from_epsg(geom.srid), CRS.from_epsg(epsg), xs, ys) + return list(zip(tx, ty)) + else: + raise ValueError("Cannot transform complex geometries to WKT") diff --git a/app/models/task.py b/app/models/task.py index e1dcd700..15779436 100644 --- a/app/models/task.py +++ b/app/models/task.py @@ -39,6 +39,7 @@ from app.cogeo import assure_cogeo from app.pointcloud_utils import is_pointcloud_georeferenced from app.testwatch import testWatch from app.security import path_traversal_check +from app.geoutils import geom_transform from nodeodm import status_codes from nodeodm.models import ProcessingNode from pyodm.exceptions import NodeResponseError, NodeConnectionError, NodeServerError, OdmError @@ -1089,6 +1090,12 @@ class Task(models.Model): } } + def get_projected_crop_bounds(self): + if self.crop is None or self.epsg is None: + return None + + return geom_transform(self.crop, self.epsg) + def get_model_display_params(self): """ Subset of a task fields used in the 3D model display view @@ -1099,7 +1106,8 @@ class Task(models.Model): 'available_assets': self.available_assets, 'public': self.public, 'public_edit': self.public_edit, - 'epsg': self.epsg + 'epsg': self.epsg, + 'crop_projected': self.get_projected_crop_bounds() } def generate_deferred_asset(self, archive, directory, stream=False): diff --git a/app/pointcloud_utils.py b/app/pointcloud_utils.py index 6777e4a7..c48fdbbf 100644 --- a/app/pointcloud_utils.py +++ b/app/pointcloud_utils.py @@ -3,7 +3,7 @@ import os import subprocess import json import rasterio -from app.api.geoutils import geom_transform_wkt_bbox +from app.geoutils import geom_transform_wkt_bbox from django.contrib.gis.geos import GEOSGeometry from app.security import double_quote diff --git a/app/raster_utils.py b/app/raster_utils.py index f93acc18..d061df0d 100644 --- a/app/raster_utils.py +++ b/app/raster_utils.py @@ -14,7 +14,7 @@ from rio_tiler.colormap import cmap as colormap, apply_cmap from rio_tiler.errors import InvalidColorMapName from app.api.hsvblend import hsv_blend from app.api.hillshade import LightSource -from app.api.geoutils import geom_transform_wkt_bbox +from app.geoutils import geom_transform_wkt_bbox from rio_tiler.io import COGReader from rasterio.warp import calculate_default_transform, reproject, Resampling diff --git a/app/static/app/js/ModelView.jsx b/app/static/app/js/ModelView.jsx index 701f3fa3..57121ad1 100644 --- a/app/static/app/js/ModelView.jsx +++ b/app/static/app/js/ModelView.jsx @@ -363,6 +363,13 @@ class ModelView extends React.Component { this.setState({error: "Could not load point cloud. This task doesn't seem to have one. Try processing the task again."}); return; } + + // Set crop vertices if needed + if (this.props.task.crop_projected && this.props.task.crop_projected.length >= 3){ + e.pointcloud.material.cropVertices = this.props.task.crop_projected.map(coord => { + return new THREE.Vector3(coord[0], coord[1], 0.0); + }); + } // Automatically load 3D model if required if (this.hasTexturedModel() && this.props.modelType === "mesh"){ diff --git a/app/static/app/js/vendor/potree/build/potree/potree.js b/app/static/app/js/vendor/potree/build/potree/potree.js index fd2ada44..db420e84 100644 --- a/app/static/app/js/vendor/potree/build/potree/potree.js +++ b/app/static/app/js/vendor/potree/build/potree/potree.js @@ -57597,6 +57597,10 @@ uniform int clipMethod; uniform mat4 uClipPolygonWVP[num_clippolygons]; #endif +#if defined(num_cropvertices) && num_cropvertices > 0 + uniform vec2 uCropVertices[num_cropvertices]; +#endif + uniform float size; uniform float minSize; @@ -58241,6 +58245,21 @@ float getPointSize(){ return pointSize; } +#if defined(num_cropvertices) && num_cropvertices > 0 +bool pointInCropPolygon(vec3 point) { + int j = num_cropvertices - 1; + bool c = false; + for (int i = 0; i < num_cropvertices; i++) { + if ( ((uCropVertices[i].y>point.y) != (uCropVertices[j].y>point.y)) && + (point.x<(uCropVertices[j].x-uCropVertices[i].x) * (point.y-uCropVertices[i].y) / (uCropVertices[j].y-uCropVertices[i].y) + uCropVertices[i].x) ){ + c = !c; + } + j = i; + } + return c; +} +#endif + #if defined(num_clippolygons) && num_clippolygons > 0 bool pointInClipPolygon(vec3 point, int polyIdx) { @@ -58353,6 +58372,16 @@ void doClipping(){ } #endif + #if defined(num_cropvertices) && num_cropvertices > 0 + // vec4 worldPos = modelMatrix * vec4(position, 1.0); + + if (!pointInCropPolygon(position)){ + gl_Position = vec4(100.0, 100.0, 100.0, 0.0); + return; + } + #endif + + #if defined(num_clippolygons) && num_clippolygons > 0 for(int i = 0; i < num_clippolygons; i++) { bool inside = pointInClipPolygon(position, i); @@ -63192,6 +63221,18 @@ void main() { } } + if (material.cropVertices && material.cropVertices.length > 0){ + let flattenedVertices = new Array(2 * material.cropVertices.length); + for(let i = 0; i < material.cropVertices.length; i++){ + // We simply subtract the translation elements from the world matrix + // since we know that the other parameters are the identity matrix + flattenedVertices[i * 2 + 0] = material.cropVertices[i].x - world.elements[12]; + flattenedVertices[i * 2 + 1] = material.cropVertices[i].y - world.elements[13]; + } + const lCropVertices = shader.uniformLocations["uCropVertices[0]"]; + gl.uniform2fv(lCropVertices, flattenedVertices); + } + //shader.setUniformMatrix4("modelMatrix", world); //shader.setUniformMatrix4("modelViewMatrix", worldView); @@ -63479,6 +63520,7 @@ void main() { let numClipBoxes = (material.clipBoxes && material.clipBoxes.length) ? material.clipBoxes.length : 0; let numClipSpheres = (params.clipSpheres && params.clipSpheres.length) ? params.clipSpheres.length : 0; let numClipPolygons = (material.clipPolygons && material.clipPolygons.length) ? material.clipPolygons.length : 0; + let numCropVertices = (material.cropVertices && material.cropVertices.length) ? material.cropVertices.length : 0; let defines = [ `#define num_shadowmaps ${shadowMaps.length}`, @@ -63486,6 +63528,7 @@ void main() { `#define num_clipboxes ${numClipBoxes}`, `#define num_clipspheres ${numClipSpheres}`, `#define num_clippolygons ${numClipPolygons}`, + `#define num_cropvertices ${numCropVertices}`, ];