Potree cropping support

app/api/geoutils.py

@@ -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")

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):

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

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")

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):

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
 

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
 

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"){

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}`,
 					];