Calculate separate isolines

[gabrielDevlog]

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);
    }
}

[mthh]

Thanks for your feedback.

Indeed, the algorithm computes closed rings.
However I think you can obtain the LineStrings you want by clipping the result by another LineString which represents the "outside" of your grid.

I didn't find how to do such an operation with the geo crate, however, maybe you can iterate on the points of your result and identify which points belong to the "outside" of the grid, an split the result on these points to reconstruct the desired output ?

See the following code, which is adapted from yours, for example (it isn't very elaborate - there must be a way to do it better or shorter)

use contour::{ContourBuilder, Line};
use std::error::Error;

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)
        .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::{Point, line_string, LineString, 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];

        let x0 = 10;
        let y0 = 10;
        let xstep = 1;
        let ystep = 1;
        let dx = 3;
        let dy = 3;

        // Calculate isolines
        let lines = isolines(x0 as f64, y0 as f64, xstep as f64, ystep as f64, dx, dy, &grid, &thresholds).unwrap();

        let multi_lines_output: MultiLineString = as_multi_line(lines)[0].clone();

        // We have a vec containing only one MultiLineString composed of only one LineString
        // so we unwrap this LineString here for simplicity :
        let result_linestring: LineString = multi_lines_output.iter().next().unwrap().clone();

        // The points that will be used to clip the returned ring
        let mut clipping_points = Vec::new();

        // Generate the possible clipping points
        for x in x0..(x0 + dx) {
            for y in y0..=(y0 + dy) {
                let _y = y as f64 - ystep as f64 / 2.0;
                clipping_points.push(
                  Point::new(x as f64, _y)
                );
            }
        }

        // We will iterate on the points of the resulting lines
        // until we encounter the clipping points (and collect the points that are
        // not the clipping point in order to reconstruct the lines)
        let mut linestrings: Vec<MultiLineString> = Vec::new();
        let mut current_points = Vec::new();

        for pt in result_linestring.points() {
            // A flag that will be true if we need to split the current linestring
            let mut found = false;
            // Iterate over the possible clipping points
            // to see if we need to split the current linestring
            for clipping_point in &clipping_points {
                if pt == *clipping_point {
                    found = true;
                    break;
                }
            }

            if !found {
                // If not found, we push the point and just continue
                current_points.push(pt);
            } else {
                // If found, we convert the existing points to a LineString
                // and we clear the vec that was storing the points
                let ls: LineString = current_points.clone().into();
                linestrings.push(ls.into());
                current_points.clear();
            }
        }

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

        assert_eq!(linestrings, expected_lines);
    }
}

Hope this helps.