Refactor various stuffs

Author: mthh

Cargo.toml

@@ -16,10 +16,20 @@ license = "MIT/Apache-2.0"
 geojson = { version = ">=0.16, <=0.24", optional = true }
 geo-types= { version = "0.7" }
 lazy_static = "1.0"
-serde_json = "^1.0"
+serde_json = { version = "^1.0", optional = true }
 rustc-hash = "1.0"
 slab = "0.4"
 
+[dev-dependencies]
+serde_json = "^1.0"
+
+[features]
+geojson = ["dep:geojson", "dep:serde_json"]
+
 [package.metadata.docs.rs]
 all-features = true
 rustdoc-args = ["--cfg", "docsrs"]
+
+[[example]]
+name = "ex"
+required-features = ["geojson"]

README.md

@@ -4,13 +4,12 @@
 [![Build status Appveyor](https://ci.appveyor.com/api/projects/status/uemh49tq7vy4uke6?svg=true)](https://ci.appveyor.com/project/mthh/contour-rs)
 [![Docs.rs version](https://docs.rs/contour/badge.svg)](https://docs.rs/contour/)
 
-Computes *isorings* and __*contour polygons*__ by applying [marching squares](https://en.wikipedia.org/wiki/Marching_squares) to a rectangular array of numeric values.  
-Outputs ring coordinates or polygons contours (represented using geo-types [MultiPolygon](https://docs.rs/geo-types/latest/geo_types/geometry/struct.MultiPolygon.html)s).
-For each threshold value, the contour polygon are representing the area where the input values are greater than or equal to the threshold value.
+Computes *isolines*, *contour polygons* and *isobands* by applying [marching squares](https://en.wikipedia.org/wiki/Marching_squares) to a rectangular array of numeric values.  
+Outputs ring coordinates or [geo-types](https://docs.rs/geo-types/) geometries.
 
-The generated contours can also easily be serialised to GeoJSON.
+The results can also easily be serialised to GeoJSON.
 
-*Note : This is a port of [d3-contour](https://github.com/d3/d3-contour).*  
+*Note : The core of the algorithm is ported from [d3-contour](https://github.com/d3/d3-contour).*  
 
 <div style="text-align:center"><a href="https://mthh.github.io/wasm_demo_contour/"><img src ="https://raw.githubusercontent.com/mthh/contour-rs/master/illustration.png" /></a></div><br>
 
@@ -31,14 +30,16 @@ extern crate contour;
 
 The API exposes:
 - a `ContourBuilder` struct, which computes isorings coordinates for a `Vec` of threshold values and transform them either :
-  - in `Contour`s (a type containing the threshold value and the geometry as a `MultiPolygon`, easily serializable to GeoJSON), or,
-  - in `Line`s (a type containing the threshold value and the geometry as a `MultiLineString`, easily serializable to GeoJSON).
-
+  - in `Contour`s (a type containing the threshold value and the geometry as a `MultiPolygon`), or,
+  - in `Line`s (a type containing the threshold value and the geometry as a `MultiLineString`).
+  - in `Band`s (a type containing a minimum value, a maximum value and the geometry as a `MultiPolygon`).
 
 - a `contour_rings` function, which computes isorings coordinates for a single threshold value (*returns a `Vec` of rings coordinates* - this is what is used internally by the `ContourBuilder`).
 
 `ContourBuilder` is the recommended way to use this crate, as it is more flexible and easier to use (it enables to specify the origin and the step of the grid, and to smooth the contours, while `contour_rings` only speak in grid coordinates and doesn't smooth the resulting rings).
 
+`Line`, `Contour` and `Band` can be serialised to GeoJSON using the `geojson` feature.
+
 ## Example:
 
 **Without defining origin and step:**
@@ -147,12 +148,13 @@ Feature {
 Demo of this crate compiled to WebAssembly and used from JavaScript : [wasm_demo_contour](https://mthh.github.io/wasm_demo_contour/).
 
 ## Difference with the [contour-isobands](https://crates.io/crates/contour-isobands) crate (from [mthh/contour-isobands-rs](https://github.com/mthh/contour-isobands-rs) repository)
-While this crate computes isolines (cf. [wikipedia:Marching_squares](https://en.wikipedia.org/wiki/Marching_squares))
-and their corresponding polygons (i.e. polygons that contain all points above the threshold defined for a given isoline),
-[contour-isobands-rs](https://github.com/mthh/contour-isobands-rs) computes isobands (cf. [wikipedia:Marching_squares#Isobands](https://en.wikipedia.org/wiki/Marching_squares#Isobands))
-and their corresponding polygons (i.e. contour polygons that contain all points between a minimum and a maximum bound).
 
-Depending on the desired use of the result, the [`contour-isobands`](https://crates.io/crates/contour-isobands) crate may be more suitable than this [`contour`](https://crates.io/crates/contour) crate.
+While this crate computes isolines (cf. [wikipedia:Marching_squares](https://en.wikipedia.org/wiki/Marching_squares)) from which it derives
+contour polygons *(i.e. polygons that contain all points above the threshold defined for a given isoline)* 
+and isobands *(i.e. polygons that contain all points between a minimum and a maximum bound)*,
+[contour-isobands-rs](https://github.com/mthh/contour-isobands-rs) is only dedicated to compute isobands and use a slightly different 
+implementation of the marching squares algorithm for the disambiguation of saddle points (cf. [wikipedia:Marching_squares](https://en.wikipedia.org/wiki/Marching_squares)).
+
 
 ## License
 

benches/bench.rs

@@ -43,13 +43,13 @@ static VALUES2: [f64; 238] = [
 ];
 
 #[bench]
-fn bench_build_geojson_contours_multiple_thresholds(b: &mut Bencher) {
+fn bench_build_contours_multiple_thresholds(b: &mut Bencher) {
     let c = ContourBuilder::new(14, 17, true);
     b.iter(|| black_box(c.contours(&VALUES2, &[0.5, 1.5, 2.5])));
 }
 
 #[bench]
-fn bench_build_geojson_contours_multiple_thresholds_and_x_y_steps_and_origins(b: &mut Bencher) {
+fn bench_build_contours_multiple_thresholds_and_x_y_steps_and_origins(b: &mut Bencher) {
     let c = ContourBuilder::new(14, 17, true)
         .x_step(0.5)
         .y_step(0.5)
@@ -79,3 +79,59 @@ fn bench_build_isoring(b: &mut Bencher) {
 fn bench_build_isoring_values2(b: &mut Bencher) {
     b.iter(|| black_box(contour_rings(&VALUES2, 1.5, 14, 17)));
 }
+
+#[bench]
+fn bench_contourbuilder_isobands_volcano_without_xy_step_xy_origin(b: &mut Bencher) {
+    let data_str = include_str!("../tests/fixtures/volcano.json");
+    let raw_data: serde_json::Value = serde_json::from_str(data_str).unwrap();
+    let matrix: Vec<f64> = raw_data["data"]
+        .as_array()
+        .unwrap()
+        .iter()
+        .map(|x| x.as_f64().unwrap())
+        .collect();
+    let h = raw_data["height"].as_u64().unwrap() as u32;
+    let w = raw_data["width"].as_u64().unwrap() as u32;
+
+    b.iter(|| {
+        black_box(
+            ContourBuilder::new(w, h, true)
+                .isobands(
+                    &matrix,
+                    &[
+                        90., 95., 100., 105., 110., 115., 120., 125., 130., 135., 140., 145., 150.,
+                        155., 160., 165., 170., 175., 180., 185., 190., 195., 200.,
+                    ],
+                )
+                .unwrap(),
+        )
+    });
+}
+
+#[bench]
+fn bench_contourbuilder_isobands_pot_pop_fr_without_xy_step_xy_origin(b: &mut Bencher) {
+    let data_str = include_str!("../tests/fixtures/pot_pop_fr.json");
+    let raw_data: serde_json::Value = serde_json::from_str(data_str).unwrap();
+    let matrix: Vec<f64> = raw_data["data"]
+        .as_array()
+        .unwrap()
+        .iter()
+        .map(|x| x.as_f64().unwrap())
+        .collect();
+    let h = raw_data["height"].as_u64().unwrap() as u32;
+    let w = raw_data["width"].as_u64().unwrap() as u32;
+
+    b.iter(|| {
+        black_box(
+            ContourBuilder::new(w, h, true)
+                .isobands(
+                    &matrix,
+                    &[
+                        0.001, 105483.25, 527416.25, 1054832.5, 2109665., 3164497.5, 4219330.,
+                        5274162.5, 6328995., 7383827.5, 8438660., 9704459., 10548326.,
+                    ],
+                )
+                .unwrap(),
+        )
+    });
+}

examples/ex.rs

@@ -0,0 +1,87 @@
+use contour::ContourBuilder;
+use geojson::{FeatureCollection, GeoJson};
+use std::fs::File;
+use std::io::{BufWriter, Write};
+
+fn main() {
+    let pot_pop_fr = include_str!("../tests/fixtures/pot_pop_fr.json");
+    let raw_data: serde_json::Value = serde_json::from_str(pot_pop_fr).unwrap();
+    let matrix: Vec<f64> = raw_data["data"]
+        .as_array()
+        .unwrap()
+        .iter()
+        .map(|x| x.as_f64().unwrap())
+        .collect();
+    let h = raw_data["height"].as_u64().unwrap() as u32;
+    let w = raw_data["width"].as_u64().unwrap() as u32;
+
+    let x_origin = -6.144721171428571;
+    let y_origin = 51.78171334283718;
+    let x_step = 0.11875873095057177;
+    let y_step = -0.08993203637245273;
+
+    let contours = ContourBuilder::new(w, h, true)
+        .x_step(x_step)
+        .y_step(y_step)
+        .x_origin(x_origin)
+        .y_origin(y_origin)
+        .isobands(
+            &matrix,
+            &[
+                0., 105483.25, 310000., 527416.25, 850000., 1054832.5, 2109665., 3164497.5,
+                4219330., 5274162.5, 6328995., 7383827.5, 8438660., 9704459., 10548326.,
+            ],
+        )
+        .unwrap();
+
+    let features = contours
+        .iter()
+        .map(|contour| contour.to_geojson())
+        .collect::<Vec<geojson::Feature>>();
+
+    let geojson_str = GeoJson::from(FeatureCollection {
+        bbox: None,
+        features,
+        foreign_members: None,
+    })
+    .to_string();
+
+    let mut file_writer = BufWriter::new(File::create("/tmp/example-output.geojson").unwrap());
+    file_writer.write(&geojson_str.as_bytes()).unwrap();
+
+    let volcano = include_str!("../tests/fixtures/volcano.json");
+    let raw_data: serde_json::Value = serde_json::from_str(volcano).unwrap();
+    let matrix: Vec<f64> = raw_data["data"]
+        .as_array()
+        .unwrap()
+        .iter()
+        .map(|x| x.as_f64().unwrap())
+        .collect();
+    let h = raw_data["height"].as_u64().unwrap() as u32;
+    let w = raw_data["width"].as_u64().unwrap() as u32;
+
+    let contours = ContourBuilder::new(w, h, true)
+        .isobands(
+            &matrix,
+            &[
+                90., 95., 100., 105., 110., 115., 120., 125., 130., 135., 140., 145., 150., 155.,
+                160., 165., 170., 175., 180., 185., 190., 195., 200.,
+            ],
+        )
+        .unwrap();
+
+    let features = contours
+        .iter()
+        .map(|contour| contour.to_geojson())
+        .collect::<Vec<geojson::Feature>>();
+
+    let geojson_str = GeoJson::from(FeatureCollection {
+        bbox: None,
+        features,
+        foreign_members: None,
+    })
+    .to_string();
+
+    let mut file_writer = BufWriter::new(File::create("/tmp/example-output2.geojson").unwrap());
+    file_writer.write(&geojson_str.as_bytes()).unwrap();
+}

src/area.rs

@@ -1,4 +1,4 @@
-use crate::contour::Pt;
+use crate::Pt;
 
 pub fn area(ring: &[Pt]) -> f64 {
     let mut i = 0;

src/band.rs

@@ -0,0 +1,72 @@
+use geo_types::MultiPolygon;
+
+/// An isoband has the geometry and min / max values of a contour ring, built by [`ContourBuilder`].
+#[derive(Debug, Clone)]
+pub struct Band {
+    pub(crate) geometry: MultiPolygon,
+    pub(crate) min_v: f64,
+    pub(crate) max_v: f64,
+}
+
+impl Band {
+    /// Borrow the [`MultiPolygon`](geo_types::MultiPolygon) geometry of this contour.
+    pub fn geometry(&self) -> &MultiPolygon {
+        &self.geometry
+    }
+
+    /// Get the owned polygons and thresholds (min and max) of this band.
+    pub fn into_inner(self) -> (MultiPolygon, f64, f64) {
+        (self.geometry, self.min_v, self.max_v)
+    }
+
+    /// Get the minimum value used to construct this band.
+    pub fn min_v(&self) -> f64 {
+        self.min_v
+    }
+
+    /// Get the maximum value used to construct this band.
+    pub fn max_v(&self) -> f64 {
+        self.max_v
+    }
+
+    #[cfg(feature = "geojson")]
+    /// Convert the band to a struct from the `geojson` crate.
+    ///
+    /// To get a string representation, call to_geojson().to_string().
+    /// ```
+    /// use contour::ContourBuilder;
+    ///
+    /// let builder = ContourBuilder::new(10, 10, false);
+    /// # #[rustfmt::skip]
+    /// let contours = builder.isobands(&[
+    /// // ...ellided for brevity
+    /// #     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 1., 2., 1., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 1., 2., 1., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 1., 2., 1., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 1., 1., 1., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
+    /// #     0., 0., 0., 0., 0., 0., 0., 0., 0., 0.
+    /// ], &[0.5, 1.5, 2.5]).unwrap();
+    ///
+    /// let geojson_string = contours[0].to_geojson().to_string();
+    ///
+    /// assert_eq!(&geojson_string[0..27], r#"{"geometry":{"coordinates":"#);
+    /// ```
+    pub fn to_geojson(&self) -> geojson::Feature {
+        let mut properties = geojson::JsonObject::with_capacity(2);
+        properties.insert("min_v".to_string(), self.min_v.into());
+        properties.insert("max_v".to_string(), self.max_v.into());
+
+        geojson::Feature {
+            bbox: None,
+            geometry: Some(geojson::Geometry::from(self.geometry())),
+            id: None,
+            properties: Some(properties),
+            foreign_members: None,
+        }
+    }
+}