mthh / contour-rs Goto Github PK

I think the area function in src/area.rs is not handling the wraparound case correctly:

pub fn area(ring: &[Pt]) -> Float {
    let mut i = 0;
    let n = ring.len() - 1;
    let mut area = ring[n - 1].y * ring[0].x - ring[n - 1].x * ring[0].y;
    while i < n {
        i += 1;
        area += ring[i - 1].y * ring[i].x - ring[i - 1].x * ring[i].y;
    }
    area
}

Say ring has 10 elements; then n == 9 and so the initialization of area is pairing up ring[0] and ring[8] (rather than ring[9]).

A correct implementation could be something like:

pub fn area(poly: &[Pt]) -> Float {
    let n = poly.len();
    let mut area = poly[n - 1].y * poly[0].x - poly[n - 1].x * poly[0].y;
    for i in 1..n {
        area += poly[i - 1].y * poly[i].x - poly[i - 1].x * poly[i].y;
    }
    area
}

(For anyone else who stumbles on this issue: I think normally you'd also want to divide the result by two, but this library is only using area to detect polygon winding and relative size of areas, so it shouldn't be necessary here.)

Every time there is a new release of geojson, the types of this crate must be updated.

This isn't a huge problem - you can just keep bumping the geojson crate (82b2373, 6398263, 520f39d) as you have, but really I don't know that there's a good reason to use the types from the geojson crate as your interchange format.

The downsides are, as mentioned above, that it triggers a downstream update cycle, but also those data structures aren't very efficient ways of representing your data types.

Would you consider a PR which leverages internal types detached from the geojson crate for contours?

If you'd like, it could retain an impl Serialize that serializes to geojson or (my preference) add a to_geojson method that explicitly converts your internal types to geojson types for those people who actually want to serialize geojson.

If I pass in a raster larger than about 40000*40000 (not sure about the exact value, but 50k*50k fails) to ContourBuilder::contours, I get an error similar to this:

thread 'main' panicked at /home/nett/.cargo/registry/src/index.crates.io-6f17d22bba15001f/contour-0.12.1/src/isoringbuilder.rs:137:23:
index out of bounds: the len is 4294836225 but the index is 18446744071562067968
note: run with `RUST_BACKTRACE=1` environment variable to display a backtrace

Increasing the width of the involved local variables to i64 seems to have fixed this issue, and I'll create a PR containing those changes once I test it against real data.

Hi,
First of all, thanks you for providing this implementation.

I'm want to calculate two isolines based on following grid, with a treshold of 1:

    // [
    //   [0.3, 1.5, 0.6],
    //   [0.6, 1.5, 0.6],
    //   [0.3, 1.5, 0.6],
    // ];

We should have two vertical lines, but I end with one circular line. Seems that the code is always trying to obtain a ring. Is there a way to obtain separates lines ?

Here is a working test, that fails:

pub fn isolines(
    x0: f64,
    y0: f64,
    x_step: f64,
    y_step: f64,
    dx: u32,
    dy: u32,
    grid: &[f64],
    thresholds: &[f64],
) -> Result<Vec<Line>, Box<dyn Error>> {
    let res = ContourBuilder::new(dx, dy, true)
        // set origin 
         .x_origin(x0 - x_step / 2.)
        .y_origin(y0 - y_step / 2.)
        .x_step(x_step)
        .y_step(y_step)
        .lines(grid, thresholds)?;

    Ok(res)
}

#[cfg(test)]
mod isolines_test {
    use contour::Line;
    use geo_types::{line_string, MultiLineString};

    use super::isolines;
   
     // transform isolines output into something we can compare with geo_types
    fn as_multi_line(lines: Vec<Line>) -> Vec<MultiLineString> {
        let multi_lines_output: Vec<MultiLineString> =
            lines.iter().map(|l| l.geometry().to_owned()).collect();
        multi_lines_output
    }

 #[test]
    fn create_two_lines_3x3_with_two_vertical_isolines() {
        // Grid data
        let grid: Vec<f64> = vec![0.3, 1.5, 0.6, 0.6, 1.5, 0.6, 0.3, 1.5, 0.6];

        // Thresholds
        let thresholds: Vec<f64> = vec![1.0];

        // Calculate isolines
        let lines = isolines(10., 10., 1., 1., 3, 3, &grid, &thresholds).unwrap();

        let multi_lines_output = as_multi_line(lines);

        // We want two lines
        let expected_lines: Vec<MultiLineString> = vec![
            line_string![
                (x: 10.583333333333334, y: 10.0),
                (x: 10.444444444444445, y: 11.0),
                (x: 10.583333333333334, y: 12.0),
            ]
            .into(),
            line_string![
                (x: 11.555555555555555, y: 12.0),
                (x: 11.555555555555555, y: 11.0),
                (x: 11.555555555555555, y: 10.0),
            ]
            .into(),
        ];

        // Fails: we have only a ring
        assert_eq!(multi_lines_output, expected_lines);
    }
}

Currently, methods return an error type, which is impossible to access from other crates, and so it is impossible to make any reasonable error handling beyond "something gone wrong".

Error struct is public, but error module is not, and it is not re-exported, so we can't name this type in a code and can't match on its kind,