Add simple gridding of LAS files

Cargo.lock

@@ -154,6 +154,12 @@ dependencies = [
  "which",
 ]
 
+[[package]]
+name = "bit_field"
+version = "0.10.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "dc827186963e592360843fb5ba4b973e145841266c1357f7180c43526f2e5b61"
+
 [[package]]
 name = "bitflags"
 version = "1.3.2"
@@ -172,6 +178,12 @@ version = "3.14.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "7f30e7476521f6f8af1a1c4c0b8cc94f0bee37d91763d0ca2665f299b6cd8aec"
 
+[[package]]
+name = "bytemuck"
+version = "1.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "374d28ec25809ee0e23827c2ab573d729e293f281dfe393500e7ad618baa61c6"
+
 [[package]]
 name = "byteorder"
 version = "1.5.0"
@@ -282,6 +294,12 @@ dependencies = [
  "cc",
 ]
 
+[[package]]
+name = "color_quant"
+version = "1.1.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3d7b894f5411737b7867f4827955924d7c254fc9f4d91a6aad6b097804b1018b"
+
 [[package]]
 name = "colorchoice"
 version = "1.0.0"
@@ -313,6 +331,37 @@ dependencies = [
  "cfg-if",
 ]
 
+[[package]]
+name = "crossbeam-deque"
+version = "0.8.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "613f8cc01fe9cf1a3eb3d7f488fd2fa8388403e97039e2f73692932e291a770d"
+dependencies = [
+ "crossbeam-epoch",
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-epoch"
+version = "0.9.18"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5b82ac4a3c2ca9c3460964f020e1402edd5753411d7737aa39c3714ad1b5420e"
+dependencies = [
+ "crossbeam-utils",
+]
+
+[[package]]
+name = "crossbeam-utils"
+version = "0.8.19"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "248e3bacc7dc6baa3b21e405ee045c3047101a49145e7e9eca583ab4c2ca5345"
+
+[[package]]
+name = "crunchy"
+version = "0.2.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7a81dae078cea95a014a339291cec439d2f232ebe854a9d672b796c6afafa9b7"
+
 [[package]]
 name = "either"
 version = "1.9.0"
@@ -344,12 +393,37 @@ dependencies = [
  "windows-sys 0.52.0",
 ]
 
+[[package]]
+name = "exr"
+version = "1.71.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "832a761f35ab3e6664babfbdc6cef35a4860e816ec3916dcfd0882954e98a8a8"
+dependencies = [
+ "bit_field",
+ "flume",
+ "half",
+ "lebe",
+ "miniz_oxide",
+ "rayon-core",
+ "smallvec",
+ "zune-inflate",
+]
+
 [[package]]
 name = "fastrand"
 version = "2.0.1"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "25cbce373ec4653f1a01a31e8a5e5ec0c622dc27ff9c4e6606eefef5cbbed4a5"
 
+[[package]]
+name = "fdeflate"
+version = "0.3.4"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4f9bfee30e4dedf0ab8b422f03af778d9612b63f502710fc500a334ebe2de645"
+dependencies = [
+ "simd-adler32",
+]
+
 [[package]]
 name = "filetime"
 version = "0.2.23"
@@ -372,6 +446,15 @@ dependencies = [
  "miniz_oxide",
 ]
 
+[[package]]
+name = "flume"
+version = "0.11.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "55ac459de2512911e4b674ce33cf20befaba382d05b62b008afc1c8b57cbf181"
+dependencies = [
+ "spin",
+]
+
 [[package]]
 name = "fnv"
 version = "1.0.7"
@@ -502,6 +585,16 @@ dependencies = [
  "serde",
 ]
 
+[[package]]
+name = "gif"
+version = "0.12.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "80792593675e051cf94a4b111980da2ba60d4a83e43e0048c5693baab3977045"
+dependencies = [
+ "color_quant",
+ "weezl",
+]
+
 [[package]]
 name = "gimli"
 version = "0.28.1"
@@ -533,6 +626,15 @@ dependencies = [
  "tracing",
 ]
 
+[[package]]
+name = "half"
+version = "2.2.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "02b4af3693f1b705df946e9fe5631932443781d0aabb423b62fcd4d73f6d2fd0"
+dependencies = [
+ "crunchy",
+]
+
 [[package]]
 name = "hashbrown"
 version = "0.14.3"
@@ -664,6 +766,24 @@ dependencies = [
  "unicode-normalization",
 ]
 
+[[package]]
+name = "image"
+version = "0.24.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "034bbe799d1909622a74d1193aa50147769440040ff36cb2baa947609b0a4e23"
+dependencies = [
+ "bytemuck",
+ "byteorder",
+ "color_quant",
+ "exr",
+ "gif",
+ "jpeg-decoder",
+ "num-traits",
+ "png",
+ "qoi",
+ "tiff",
+]
+
 [[package]]
 name = "indexmap"
 version = "2.1.0"
@@ -686,6 +806,15 @@ version = "1.0.10"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "b1a46d1a171d865aa5f83f92695765caa047a9b4cbae2cbf37dbd613a793fd4c"
 
+[[package]]
+name = "jpeg-decoder"
+version = "0.3.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "f5d4a7da358eff58addd2877a45865158f0d78c911d43a5784ceb7bbf52833b0"
+dependencies = [
+ "rayon",
+]
+
 [[package]]
 name = "js-sys"
 version = "0.3.67"
@@ -732,6 +861,12 @@ version = "1.3.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "830d08ce1d1d941e6b30645f1a0eb5643013d835ce3779a5fc208261dbe10f55"
 
+[[package]]
+name = "lebe"
+version = "0.5.2"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "03087c2bad5e1034e8cace5926dec053fb3790248370865f5117a7d0213354c8"
+
 [[package]]
 name = "libc"
 version = "0.2.152"
@@ -760,6 +895,16 @@ version = "0.4.13"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "01cda141df6706de531b6c46c3a33ecca755538219bd484262fa09410c13539c"
 
+[[package]]
+name = "lock_api"
+version = "0.4.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3c168f8615b12bc01f9c17e2eb0cc07dcae1940121185446edc3744920e8ef45"
+dependencies = [
+ "autocfg",
+ "scopeguard",
+]
+
 [[package]]
 name = "log"
 version = "0.4.20"
@@ -791,6 +936,7 @@ source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "e7810e0be55b428ada41041c41f32c9f1a42817901b4ccf45fa3d4b6561e74c7"
 dependencies = [
  "adler",
+ "simd-adler32",
 ]
 
 [[package]]
@@ -941,6 +1087,19 @@ version = "0.3.29"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "2900ede94e305130c13ddd391e0ab7cbaeb783945ae07a279c268cb05109c6cb"
 
+[[package]]
+name = "png"
+version = "0.17.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1f6c3c3e617595665b8ea2ff95a86066be38fb121ff920a9c0eb282abcd1da5a"
+dependencies = [
+ "bitflags 1.3.2",
+ "crc32fast",
+ "fdeflate",
+ "flate2",
+ "miniz_oxide",
+]
+
 [[package]]
 name = "prettyplease"
 version = "0.2.16"
@@ -995,6 +1154,7 @@ dependencies = [
  "clap",
  "futures",
  "geo-types",
+ "image",
  "las",
  "proj",
  "reqwest",
@@ -1003,6 +1163,15 @@ dependencies = [
  "tokio-util",
 ]
 
+[[package]]
+name = "qoi"
+version = "0.4.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "7f6d64c71eb498fe9eae14ce4ec935c555749aef511cca85b5568910d6e48001"
+dependencies = [
+ "bytemuck",
+]
+
 [[package]]
 name = "quote"
 version = "1.0.35"
@@ -1012,6 +1181,26 @@ dependencies = [
  "proc-macro2",
 ]
 
+[[package]]
+name = "rayon"
+version = "1.8.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "fa7237101a77a10773db45d62004a272517633fbcc3df19d96455ede1122e051"
+dependencies = [
+ "either",
+ "rayon-core",
+]
+
+[[package]]
+name = "rayon-core"
+version = "1.12.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1465873a3dfdaa8ae7cb14b4383657caab0b3e8a0aa9ae8e04b044854c8dfce2"
+dependencies = [
+ "crossbeam-deque",
+ "crossbeam-utils",
+]
+
 [[package]]
 name = "redox_syscall"
 version = "0.4.1"
@@ -1130,6 +1319,12 @@ dependencies = [
  "windows-sys 0.52.0",
 ]
 
+[[package]]
+name = "scopeguard"
+version = "1.2.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "94143f37725109f92c262ed2cf5e59bce7498c01bcc1502d7b9afe439a4e9f49"
+
 [[package]]
 name = "security-framework"
 version = "2.9.2"
@@ -1202,6 +1397,12 @@ version = "1.2.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "a7cee0529a6d40f580e7a5e6c495c8fbfe21b7b52795ed4bb5e62cdf92bc6380"
 
+[[package]]
+name = "simd-adler32"
+version = "0.3.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "d66dc143e6b11c1eddc06d5c423cfc97062865baf299914ab64caa38182078fe"
+
 [[package]]
 name = "slab"
 version = "0.4.9"
@@ -1211,6 +1412,12 @@ dependencies = [
  "autocfg",
 ]
 
+[[package]]
+name = "smallvec"
+version = "1.13.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "e6ecd384b10a64542d77071bd64bd7b231f4ed5940fba55e98c3de13824cf3d7"
+
 [[package]]
 name = "socket2"
 version = "0.5.5"
@@ -1221,6 +1428,15 @@ dependencies = [
  "windows-sys 0.48.0",
 ]
 
+[[package]]
+name = "spin"
+version = "0.9.8"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6980e8d7511241f8acf4aebddbb1ff938df5eebe98691418c4468d0b72a96a67"
+dependencies = [
+ "lock_api",
+]
+
 [[package]]
 name = "strsim"
 version = "0.10.0"
@@ -1303,6 +1519,17 @@ dependencies = [
  "syn",
 ]
 
+[[package]]
+name = "tiff"
+version = "0.9.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ba1310fcea54c6a9a4fd1aad794ecc02c31682f6bfbecdf460bf19533eed1e3e"
+dependencies = [
+ "flate2",
+ "jpeg-decoder",
+ "weezl",
+]
+
 [[package]]
 name = "tinyvec"
 version = "1.6.0"
@@ -1555,6 +1782,12 @@ dependencies = [
  "wasm-bindgen",
 ]
 
+[[package]]
+name = "weezl"
+version = "0.1.7"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9193164d4de03a926d909d3bc7c30543cecb35400c02114792c2cae20d5e2dbb"
+
 [[package]]
 name = "which"
 version = "4.4.2"
@@ -1728,3 +1961,12 @@ dependencies = [
  "linux-raw-sys",
  "rustix",
 ]
+
+[[package]]
+name = "zune-inflate"
+version = "0.2.54"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "73ab332fe2f6680068f3582b16a24f90ad7096d5d39b974d1c0aff0125116f02"
+dependencies = [
+ "simd-adler32",
+]

Cargo.toml

@@ -17,3 +17,4 @@ anyhow = "1.0"
 futures = "0.3"
 bytes = "1.5"
 las = {version="0.8", features=["laz"]}
+image = "0.24"

src/main.rs

@@ -1,4 +1,15 @@
-use {clap::Parser, geo_types::Point, las::Read as LasRead, proj::Proj, tokio::io::AsyncReadExt};
+use {
+    clap::Parser,
+    geo_types::Point,
+    image::{ImageBuffer, Luma},
+    las::Read as LasRead,
+    proj::Proj,
+    std::{
+        ops::{Add, Div, Mul, Sub},
+        path::PathBuf,
+    },
+    tokio::io::AsyncReadExt,
+};
 
 mod geonb;
 
@@ -13,9 +24,207 @@ struct Args {
     #[clap(long)]
     longitude: Option<f64>,
 
-    // Print extra debug info when initializing Vulkan
+    // Print extra debug info when initializing raphics
     #[clap(long)]
     debug_init: bool,
+
+    #[clap(long)]
+    laz_file: Option<PathBuf>,
+
+    #[clap(long)]
+    grid_cell_size: Option<f64>,
+}
+
+trait FromFloatFloor<F> {
+    fn from_float_floor(f: F) -> Self;
+}
+
+trait ToFloat<F> {
+    fn to_float(self) -> F;
+}
+
+trait Integer:
+    Copy
+    + Default
+    + Add<Self, Output = Self>
+    + Sub<Self, Output = Self>
+    + Mul<Self, Output = Self>
+    + Div<Self, Output = Self>
+    + FromFloatFloor<f32>
+    + FromFloatFloor<f64>
+    + ToFloat<f32>
+    + ToFloat<f64>
+{
+    const MIN: Self;
+    const MAX: Self;
+}
+
+macro_rules! impl_int {
+    ( $t:ident ) => {
+        impl FromFloatFloor<f32> for $t {
+            fn from_float_floor(f: f32) -> Self {
+                f as $t
+            }
+        }
+        impl FromFloatFloor<f64> for $t {
+            fn from_float_floor(f: f64) -> Self {
+                f as $t
+            }
+        }
+        impl ToFloat<f32> for $t {
+            fn to_float(self) -> f32 {
+                self as f32
+            }
+        }
+        impl ToFloat<f64> for $t {
+            fn to_float(self) -> f64 {
+                self as f64
+            }
+        }
+        impl Integer for $t {
+            const MIN: $t = $t::MIN;
+            const MAX: $t = $t::MAX;
+        }
+    };
+}
+
+impl_int!(u8);
+impl_int!(i8);
+impl_int!(u16);
+impl_int!(i16);
+impl_int!(u32);
+impl_int!(i32);
+impl_int!(u64);
+impl_int!(i64);
+
+trait Unsigned {}
+
+impl Unsigned for u8 {}
+impl Unsigned for u16 {}
+impl Unsigned for u32 {}
+impl Unsigned for u64 {}
+
+trait Float:
+    Copy
+    + Default
+    + Add<Self, Output = Self>
+    + Sub<Self, Output = Self>
+    + Mul<Self, Output = Self>
+    + Div<Self, Output = Self>
+{
+    fn to_int_floor<I: Integer>(self) -> I;
+}
+
+macro_rules! impl_float {
+    ( $t:ident ) => {
+        impl Float for $t {
+            fn to_int_floor<I: Integer>(self) -> I {
+                I::from_float_floor(self)
+            }
+        }
+    };
+}
+
+impl_float!(f32);
+impl_float!(f64);
+
+struct Grid<T> {
+    width: usize,
+    height: usize,
+    data: Vec<T>,
+}
+
+impl<T> Grid<T> {
+    fn calculate_index(&self, i: usize, j: usize) -> usize {
+        assert!(i < self.width);
+        assert!(j < self.height);
+        i + j * self.width
+    }
+}
+
+impl<T: Default + Copy> Grid<T> {
+    fn new(width: usize, height: usize) -> Self {
+        Self {
+            width,
+            height,
+            data: vec![T::default(); width * height],
+        }
+    }
+}
+
+impl<T: Copy> Grid<T> {
+    fn set_cell_value(&mut self, i: usize, j: usize, value: T) {
+        let index = self.calculate_index(i, j);
+        self.data[index] = value;
+    }
+
+    fn get_cell_value(&self, i: usize, j: usize) -> T {
+        let index = self.calculate_index(i, j);
+        self.data[index]
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn get_cell_value_returns_zero_after_new() {
+        let width = 11;
+        let height = 13;
+        let target: Grid<i32> = Grid::new(width, height);
+        for i in 0..width {
+            for j in 0..height {
+                assert_eq!(0, target.get_cell_value(i, j));
+            }
+        }
+    }
+
+    #[test]
+    fn get_cell_value_returns_values_set_by_set_cell_value() {
+        fn unique_value_for_cell(i: usize, j: usize) -> i32 {
+            let high_bits = i << 12;
+            let low_bits = j;
+            assert_eq!(0, high_bits & low_bits);
+            let bits = low_bits | high_bits;
+            assert_eq!(bits, bits & 0x7fffffff);
+            bits as i32
+        }
+        let width = 11;
+        let height = 13;
+        let mut target: Grid<i32> = Grid::new(width, height);
+        for i in 0..width {
+            for j in 0..height {
+                target.set_cell_value(i, j, unique_value_for_cell(i, j));
+                assert_eq!(unique_value_for_cell(i, j), target.get_cell_value(i, j));
+            }
+        }
+        for i in 0..width {
+            for j in 0..height {
+                assert_eq!(unique_value_for_cell(i, j), target.get_cell_value(i, j));
+            }
+        }
+    }
+}
+
+fn integerize<F, I>(grid_in: Grid<F>, min: F, max: F) -> Grid<I>
+where
+    I: Integer + ToFloat<F>,
+    F: Float,
+{
+    let mut grid_out = Grid::new(grid_in.width, grid_in.height);
+    let range = max - min;
+    for i in 0..(grid_in.width as usize) {
+        for j in 0..(grid_in.height as usize) {
+            grid_out.set_cell_value(
+                i,
+                j,
+                (((grid_in.get_cell_value(i, j) - min) / range) * (I::MAX).to_float())
+                    .to_int_floor(),
+            );
+        }
+    }
+    grid_out
 }
 
 #[tokio::main]
@@ -50,5 +259,37 @@ async fn main() -> Result<(), anyhow::Error> {
         }
     }
 
+    if let Some(laz_filename) = args.laz_file {
+        let mut las_reader = las::Reader::from_path(laz_filename).unwrap();
+        let bounds = las_reader.header().bounds();
+        if let Some(grid_cell_size) = args.grid_cell_size {
+            let num_cells_x = ((bounds.max.x - bounds.min.x) / grid_cell_size).ceil() as usize;
+            let num_cells_y = ((bounds.max.y - bounds.min.y) / grid_cell_size).ceil() as usize;
+            let mut height_grid = Grid::new(num_cells_y, num_cells_x);
+            let mut weight_grid = Grid::new(num_cells_y, num_cells_x);
+            println!("Creating {} by {} grid.", num_cells_y, num_cells_x);
+
+            for wrapped_point in las_reader.points() {
+                let point = wrapped_point.unwrap();
+                let bin_i = ((point.x - bounds.min.x) / grid_cell_size).floor() as usize;
+                let bin_j = ((point.y - bounds.min.y) / grid_cell_size).floor() as usize;
+                let height: f64 = height_grid.get_cell_value(bin_i, bin_j);
+                let weight: f64 = weight_grid.get_cell_value(bin_i, bin_j);
+                height_grid.set_cell_value(
+                    bin_i,
+                    bin_j,
+                    (height * weight + point.z) / (weight + 1.0),
+                );
+                weight_grid.set_cell_value(bin_i, bin_j, weight + 1.0);
+            }
+
+            let grid_u8: Grid<u8> = integerize(height_grid, bounds.min.z, bounds.max.z);
+            let image: ImageBuffer<Luma<u8>, _> =
+                ImageBuffer::from_raw(num_cells_x as u32, num_cells_y as u32, grid_u8.data)
+                    .unwrap();
+            image.save("test.png")?;
+        }
+    }
+
     Ok(())
 }