Add a native mesh rasterization (#36)

* Add native mesh rasterization

Author: pomadchin

CHANGELOG.md

@@ -5,6 +5,11 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.0.0/),
 and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
 
 ## [Unreleased]
+
+## [2.1.5-RC1] - 2020-03-09
+### Added
+- Add a native mesh rasterization [#36](https://github.com/PDAL/java/pull/36)
+
 ### Changed
 - Move Scala syntax into the pipeline.syntax package [#33](https://github.com/PDAL/java/issues/33)
 
@@ -88,7 +93,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Changed
 - Moved from the PDAL repo and established own lifecycle.
 
-[Unreleased]: https://github.com/PDAL/java/compare/2.1.4...HEAD
+[Unreleased]: https://github.com/PDAL/java/compare/2.1.5-RC1...HEAD
+[2.1.5-RC1]: https://github.com/PDAL/java/compare/2.1.4...2.1.5-RC1
 [2.1.4]: https://github.com/PDAL/java/compare/2.1.3...2.1.4
 [2.1.3]: https://github.com/PDAL/java/compare/2.1.2...2.1.3
 [2.1.2]: https://github.com/PDAL/java/compare/2.0.0...2.1.2

core/src/main/scala/io/pdal/PointView.scala

@@ -145,13 +145,23 @@ class PointView extends Native {
   def getY(packedPoint: Array[Byte]): Double = getDouble(packedPoint, DimType.Y)
   def getZ(packedPoint: Array[Byte]): Double = getDouble(packedPoint, DimType.Z)
 
+  def getTriangularMesh(): TriangularMesh = getTriangularMesh("")
+
+  def rasterizeTriangularMesh(extent: Array[Double], cols: Int, rows: Int): Array[Double] =
+    rasterizeTriangularMesh(extent, cols, rows, DimType.Z)
+  def rasterizeTriangularMesh(extent: Array[Double], cols: Int, rows: Int, dim: String): Array[Double] =
+    rasterizeTriangularMesh(extent, cols, rows, findDimType(dim))
+  def rasterizeTriangularMesh(extent: Array[Double], cols: Int, rows: Int, dim: DimType): Array[Double] =
+    rasterizeTriangularMesh(extent, cols, rows, dim, "")
+
   @native def layout(): PointLayout
   @native def size(): Int
   @native def empty(): Boolean
   @native def getCrsProj4(): String
   @native def getCrsWKT(pretty: Boolean): String
   @native def getPackedPoint(idx: Long, dims: Array[DimType]): Array[Byte]
   @native def getPackedPoints(dims: Array[DimType]): Array[Byte]
-  @native def getTriangularMesh(name: String = ""): TriangularMesh
+  @native def getTriangularMesh(name: String): TriangularMesh
+  @native def rasterizeTriangularMesh(extent: Array[Double], cols: Int, rows: Int, dim: DimType, name: String): Array[Double]
   @native def close(): Unit
 }

core/src/test/scala/io/pdal/PipelineSpec.scala

@@ -242,5 +242,33 @@ class PipelineSpec extends TestEnvironmentSpec {
       pv.close()
       pvi.close()
     }
+
+    it("should rasterize mesh") {
+      val pvi = pipelineDelaunay.getPointViews()
+      val pv = pvi.next()
+
+      val raster = pv.rasterizeTriangularMesh(Array(635619.85, 848899.7, 638982.55, 853535.43), 100, 100)
+
+      val (mi, ma) = {
+        val sorted = raster.filter(!_.isNaN).sorted
+        sorted.head -> sorted.last
+      }
+
+      mi shouldBe 406.915 +- 1e-3
+      ma shouldBe 582.022 +- 1e-3
+
+      pv.close()
+      pvi.close()
+    }
+
+    it("should not rasterize mesh if it was not generated") {
+      val pvi = pipeline.getPointViews()
+      val pv = pvi.next()
+
+      intercept[ExecutionException] { pv.rasterizeTriangularMesh(Array(), 0, 0) }
+
+      pv.close()
+      pvi.close()
+    }
   }
 }

native/src/include/io_pdal_PointView.h

@@ -33,6 +33,7 @@
 
 #include <stdio.h>
 #include <vector>
+#include <float.h>
 #include "io_pdal_PointView.h"
 #include "JavaPipeline.hpp"
 #include "JavaTriangularMeshIterator.hpp"
@@ -52,6 +53,8 @@ using pdal::Dimension::Id;
 using pdal::PointId;
 using pdal::DimTypeList;
 using pdal::SpatialReference;
+using pdal::TriangularMesh;
+using pdal::Triangle;
 using pdal::DimType;
 using pdal::pdal_error;
 
@@ -225,6 +228,172 @@ JNIEXPORT jobject JNICALL Java_io_pdal_PointView_getTriangularMesh
     return mi;
 }
 
+JNIEXPORT jdoubleArray JNICALL Java_io_pdal_PointView_rasterizeTriangularMesh
+  (JNIEnv *env, jobject obj, jdoubleArray extent, jint cols, jint rows, jobject jDimType, jstring name)
+{
+    std::string cname = std::string(env->GetStringUTFChars(name, 0));
+
+    PointViewRawPtr *pvrp = getHandle<PointViewRawPtr>(env, obj);
+    PointViewPtr pv = pvrp->shared_pointer;
+    TriangularMesh *mesh = pv->mesh(cname);
+
+    if(mesh == NULL)
+    {
+        throwExecutionException(env, "No mesh was generated. Check that the appropriate filter is a part of a PDAL Pipeline.");
+        return env->NewDoubleArray(0);
+    }
+
+    jclass cDimType = env->GetObjectClass(jDimType);
+    jfieldID fid = env->GetFieldID(cDimType, "id", "Ljava/lang/String;");
+    jstring jid = reinterpret_cast<jstring>(env->GetObjectField(jDimType, fid));
+    
+    PointLayoutPtr pl = pv->layout();
+    DimType dimType = pl->findDimType(std::string(env->GetStringUTFChars(jid, 0)));
+    Id dimId = dimType.m_id;
+
+    int size = mesh->size();
+
+    int length = cols * rows;
+
+    jdouble *ext = env->GetDoubleArrayElements(extent, 0);
+    
+    double exmin = ext[0];
+    double eymin = ext[1];
+    double exmax = ext[2];
+    double eymax = ext[3];
+
+    env->ReleaseDoubleArrayElements(extent, ext, JNI_ABORT);
+
+    double w = (exmax - exmin) / cols;
+    double h = (eymax - eymin) / rows;
+
+    jdoubleArray result = env->NewDoubleArray(length);
+
+    double buffer[length];
+    std::fill(buffer, buffer + sizeof(buffer) / sizeof(double), strtod("NaN", NULL));
+
+    for (int id = 0; id < size; id++)
+    {
+        const Triangle& tri = (*mesh)[id];
+        PointId a = tri.m_a;
+        PointId b = tri.m_b;
+        PointId c = tri.m_c;
+
+        double v1x = pv->getFieldAs<double>(Id::X, a);
+        double v1y = pv->getFieldAs<double>(Id::Y, a);
+        double v1z = pv->getFieldAs<double>(dimId, a);
+        double s1x = pv->getFieldAs<double>(Id::X, c);
+        double s1y = pv->getFieldAs<double>(Id::Y, c);
+
+        double v2x = pv->getFieldAs<double>(Id::X, b);
+        double v2y = pv->getFieldAs<double>(Id::Y, b);
+        double v2z = pv->getFieldAs<double>(dimId, b);
+        double s2x = pv->getFieldAs<double>(Id::X, a);
+        double s2y = pv->getFieldAs<double>(Id::Y, a);
+
+        double v3x = pv->getFieldAs<double>(Id::X, c);
+        double v3y = pv->getFieldAs<double>(Id::Y, c);
+        double v3z = pv->getFieldAs<double>(dimId, c);
+        double s3x = pv->getFieldAs<double>(Id::X, b);
+        double s3y = pv->getFieldAs<double>(Id::Y, b);
+
+        double determinant = (v2y - v3y) * (v1x - v3x) + (v3x - v2x) * (v1y - v3y);
+
+        double ymin = std::min(v1y, std::min(v2y, v3y));
+        double ymax = std::max(v1y, std::max(v2y, v3y));
+
+        double scanrow0 = std::max(std::ceil((ymin - eymin) / h - 0.5), 0.0);
+        double scany = eymin + scanrow0 * h + h / 2.0;
+
+        while (scany < eymax && scany < ymax) 
+        {
+            // get x at y for edge
+            double xmin = -DBL_MAX;
+            double xmax = DBL_MAX;
+
+            if(s1y != v1y) 
+            {
+                double t = (scany - v1y) / (s1y - v1y);
+                double xAtY1 = v1x + t * (s1x - v1x);
+
+                if(v1y < s1y) 
+                {
+                    // Lefty
+                    if(xmin < xAtY1) { xmin = xAtY1; }
+                } 
+                else 
+                {
+                    // Righty
+                    if(xAtY1 < xmax) { xmax = xAtY1; }
+                }
+            }
+
+            if(s2y != v2y) 
+            {
+                double t = (scany - v2y) / (s2y - v2y);
+                double xAtY2 = v2x + t * (s2x - v2x);
+
+                if(v2y < s2y) 
+                {
+                    // Lefty
+                    if(xmin < xAtY2) { xmin = xAtY2; }
+                } 
+                else 
+                {
+                    // Righty
+                    if(xAtY2 < xmax) { xmax = xAtY2; }
+                }
+            }
+
+            if(s3y != v3y) 
+            {
+                double t = (scany - v3y) / (s3y - v3y);
+                double xAtY3 = v3x + t * (s3x - v3x);
+
+                if(v3y < s3y) 
+                {
+                    // Lefty
+                    if(xmin < xAtY3) { xmin = xAtY3; }
+                } 
+                else 
+                {
+                    // Righty
+                    if(xAtY3 < xmax) { xmax = xAtY3; }
+                }
+            }
+
+            double scancol0 = std::max(std::ceil((xmin - exmin) / w - 0.5), 0.0);
+            double scanx = exmin + scancol0 * w + w / 2;
+
+            while (scanx < exmax && scanx < xmax) 
+            {
+                int col = (int)((scanx - exmin) / w);
+                int row = (int)((eymax - scany) / h);
+
+                if(0 <= col && col < cols && 0 <= row && row < rows)
+                {
+
+                    double lambda1 = ((v2y - v3y) * (scanx - v3x) + (v3x - v2x) * (scany - v3y)) / determinant;
+                    double lambda2 = ((v3y - v1y) * (scanx - v3x) + (v1x - v3x) * (scany - v3y)) / determinant;
+                    double lambda3 = 1.0 - lambda1 - lambda2;
+
+                    double z = lambda1 * v1z + lambda2 * v2z + lambda3 * v3z;
+
+                    buffer[row * cols + col] = z;
+                }
+
+                scanx += w;
+            }
+
+            scany += h;
+        }
+    }
+
+    env->SetDoubleArrayRegion(result, 0, length, buffer);
+
+    return result;
+}
+
 JNIEXPORT void JNICALL Java_io_pdal_PointView_close
   (JNIEnv *env, jobject obj)
 {