x_origin and y_origin are off by 1/2 a step size

[pfbuxton]

Hello. This is a nice and useful crate.

I found that x_origin and y_origin are off by 1/2 a step size. You could either change how the API works or keep it as is and add extra documentation.

Here is a reproducible example, the important parts are .x_origin(x[0] - d_x/2.0) and .y_origin(y[0] - d_y/2.0)

src/main.rs

use contour::ContourBuilder;
use ndarray::{Array1, Array2};


fn main() {
    let n_x: usize = 10;
    let n_y: usize = 9;
    let x: Array1<f64> = ndarray::Array1::linspace(0.0, 1.0, n_x);
    let y: Array1<f64> = ndarray::Array1::linspace(0.0, 1.0, n_y);
    let d_x = x[1] - x[0];
    let d_y = y[1] - y[0];

    let mut h: Array2<f64> = Array2::zeros((n_y, n_x));

    for i_x in 0..n_x {
        for i_y in 0..n_y {
            h[[i_y, i_x]] = x[i_x] * y[i_y];
        }
    }

    let c: ContourBuilder = ContourBuilder::new(n_x, n_y, true) // x dim., y dim., smoothing
        .x_step(d_x)
        .y_step(d_y)
        .x_origin(x[0] - d_x/2.0)
        .y_origin(y[0] - d_y/2.0);

    let h_flattened: Vec<f64> = h.iter().cloned().collect();

    let contour_results: Vec<contour::Contour> = c.contours(&h_flattened, &[0.5]).unwrap(); // values, thresholds

    let first_contour = contour_results[0].geometry(); // The [0] is because we have only supplied one threshold

    // Start with an empty Vec
    let mut x_contour: Vec<f64> = Vec::new();
    let mut y_contour: Vec<f64> = Vec::new();

    for polygon in first_contour.iter() {
        for coord in polygon.exterior().coords() {
            x_contour.push(coord.x);
            y_contour.push(coord.y);
        }
    }

    let x_contour: Array1<f64> = Array1::from_vec(x_contour);
    let y_contour: Array1<f64> = Array1::from_vec(y_contour);

    println!("{}", x_contour);
    println!("{}", y_contour);

}

Cargo.toml

[package]
name = "contour_example"
version = "0.1.0"
edition = "2021"

[dependencies]
contour = "0.13.1"
ndarray = "0.16.1"

and comparing with Python (doing a copy and paste from the rust prtinln):

import matplotlib.pyplot as plt
import numpy as np

n_x = 10
n_y = 9
x = np.linspace(0.0, 1.0, n_x)
y = np.linspace(0.0, 1.0, n_y)

h = np.zeros([n_y, n_x])
for i_x in range(0, n_x):
    for i_y in range(0, n_y):
        h[i_y, i_x] = x[i_x] * y[i_y]

plt.contour(x, y, h, [0.5], colors="blue")

# Plot the rust result
x_rust = [1.0555555555555556, 1.0555555555555556, 1.0555555555555556, 1.0555555555555556, 1.0555555555555556, 1, 1, 0.8888888888888888, 0.8, 0.7777777777777777, 0.6666666666666666, 0.6666666666666666, 0.5714285714285713, 0.5555555555555555, 0.5, 0.5555555555555555, 0.6666666666666666, 0.7777777777777777, 0.8888888888888888, 1, 1.0555555555555556]
y_rust = [1, 0.875, 0.75, 0.625, 0.5, 0.5, 0.5, 0.5625, 0.625, 0.6428571428571429, 0.75, 0.75, 0.875, 0.8999999999999999, 1, 1.0625, 1.0625, 1.0625, 1.0625, 1.0625, 1]
plt.plot(x_rust, y_rust, linestyle="--", color="red")

plt.savefig('python_vs_rust.svg')

and we get this figure:

There is something a bit odd with the points which are at the boundary as they are now at the boundary + 1/2 step size, which is now outside the grid!

Thanks, Peter