inpoly/
lib.rs

1use ndarray::Array1;
2use ndarray::Array2;
3use ndarray::Axis;
4use ndarray_stats::QuantileExt;
5
6fn searchsorted_left(y_column: &[f64], y_values: &[f64]) -> Vec<usize> {
7    y_values
8        .iter()
9        .map(|&y| {
10            y_column
11                .binary_search_by(|&val| val.partial_cmp(&y).unwrap())
12                .unwrap_or_else(|x| x)
13        })
14        .collect()
15}
16
17fn searchsorted_right(y_column: &[f64], y_values: &[f64]) -> Vec<usize> {
18    y_values
19        .iter()
20        .map(|&y| {
21            match y_column.binary_search_by(|&val| val.partial_cmp(&y).unwrap()) {
22                Ok(pos) => {
23                    // Move to the next position for right side
24                    let mut pos = pos;
25                    while pos < y_column.len() && y_column[pos] == y {
26                        pos += 1;
27                    }
28                    pos
29                }
30                Err(pos) => pos,
31            }
32        })
33        .collect()
34}
35
36#[allow(non_snake_case)]
37fn _inpoly(
38    vert: &Array2<f64>,
39    node: &Array2<f64>,
40    edge: &Array2<usize>,
41    ftol: &f64,
42    lbar: &f64,
43) -> (Array1<bool>, Array1<bool>) {
44    // verbatim from the Python code, not entirely sure why the use of powf here.
45    let feps = ftol * (lbar.powf(1.0));
46    let veps = ftol * (lbar.powf(1.0));
47
48    let vnum = vert.nrows();
49    let mut stat = Array1::from_elem(vnum, false);
50    let mut bnds = Array1::from_elem(vnum, false);
51
52    // compute y-range overlap
53    let XONE = node
54        .select(Axis(0), &edge.column(0).to_owned().into_raw_vec())
55        .column(0)
56        .to_owned();
57    let XTWO = node
58        .select(Axis(0), &edge.column(1).to_owned().into_raw_vec())
59        .column(0)
60        .to_owned();
61    let YONE = node
62        .select(Axis(0), &edge.column(0).to_owned().into_raw_vec())
63        .column(1)
64        .to_owned();
65    let YTWO = node
66        .select(Axis(0), &edge.column(1).to_owned().into_raw_vec())
67        .column(1)
68        .to_owned();
69    let XMIN: Vec<f64> = XONE
70        .iter()
71        .zip(XTWO.iter())
72        .map(|(a, b)| a.min(*b))
73        .collect();
74    let XMAX: Vec<f64> = XONE
75        .iter()
76        .zip(XTWO.iter())
77        .map(|(a, b)| a.max(*b))
78        .collect();
79
80    let XMIN: Vec<f64> = XMIN.iter().map(|x| x - veps).collect();
81    let XMAX: Vec<f64> = XMAX.iter().map(|x| x + veps).collect();
82    let YMIN = &YONE.mapv(|y| y - veps);
83    let YMAX = &YTWO.mapv(|y| y + veps);
84
85    let YDEL = &YTWO - &YONE;
86    let XDEL = &XTWO - &XONE;
87
88    let EDEL = &XDEL.mapv(|x| x.abs()) + &YDEL;
89
90    // Get the sorted indices based on the second column (y values)
91    let y_column = vert.column(1).to_vec();
92    let mut ivec: Vec<usize> = (0..vert.nrows()).collect();
93    ivec.sort_by(|&a, &b| y_column[a].partial_cmp(&y_column[b]).unwrap());
94    let ione = searchsorted_left(&y_column, &YMIN.to_vec());
95    let itwo = searchsorted_right(&y_column, &YMAX.to_vec());
96
97    for epos in 0..edge.nrows() {
98        let xone = XONE[epos];
99        let xtwo = XTWO[epos];
100        let yone = YONE[epos];
101        let ytwo = YTWO[epos];
102        let xmin = XMIN[epos];
103        let xmax = XMAX[epos];
104        let edel = EDEL[epos];
105        let xdel = XDEL[epos];
106        let ydel = YDEL[epos];
107
108        for jpos in ione[epos]..itwo[epos] {
109            let jvrt = ivec[jpos];
110            if bnds[jvrt] == true {
111                continue;
112            }
113            let xpos = vert[[jvrt, 0]];
114            let ypos = vert[[jvrt, 1]];
115            if xpos >= xmin {
116                if xpos <= xmax {
117                    // compute crossing number
118                    let mul1 = ydel * (xpos - xone);
119                    let mul2 = xdel * (ypos - yone);
120                    if (feps * edel) >= (mul2 - mul1).abs() {
121                        // BNDS -- approx. on edge
122                        bnds[jvrt] = true;
123                        stat[jvrt] = true;
124                    } else if (ypos == yone) && (xpos == xone) {
125                        // BNDS -- match about ONE
126                        bnds[jvrt] = true;
127                        stat[jvrt] = true;
128                    } else if (ypos == ytwo) && (xpos == xtwo) {
129                        // BNDS -- match about TWO
130                        bnds[jvrt] = true;
131                        stat[jvrt] = true;
132                    } else if (mul1 <= mul2) && (ypos >= yone) && (ypos < ytwo) {
133                        stat[jvrt] = !stat[jvrt]
134                    }
135                };
136            } else if (ypos >= yone) && (ypos < ytwo) {
137                // advance crossing number
138                stat[jvrt] = !stat[jvrt]
139            }
140        }
141    }
142    (stat, bnds)
143}
144
145#[allow(non_snake_case)]
146pub fn inpoly2(
147    vert: &Array2<f64>,
148    node: &Array2<f64>,
149    edge: Option<&Array2<usize>>,
150    ftol: Option<&f64>,
151) -> (Array1<bool>, Array1<bool>) {
152    let vnum = vert.nrows();
153
154    let mut STAT = Array1::from_elem(vnum, false);
155    let mut BNDS = Array1::from_elem(vnum, false);
156
157    if node.is_empty() {
158        return (STAT, BNDS);
159    }
160    let edge = match edge {
161        None => {
162            let n = node.nrows();
163            let mut e = Array2::zeros((n, 2));
164            for i in 0..n - 1 {
165                e[[i, 0]] = i as usize;
166                e[[i, 1]] = (i + 1) as usize;
167            }
168            e[[n - 1, 0]] = (n - 1) as usize;
169            Some(e)
170        }
171        Some(edge) => Some(edge.clone()),
172    };
173
174    let edge = edge.unwrap();
175    let used = edge.iter().map(|&e| e as usize).collect::<Vec<_>>();
176    let sel0 = node.select(Axis(0), &used);
177    let col0 = sel0.column(0);
178    let xmin = col0.min().unwrap();
179    let xmax = col0.max().unwrap();
180    let col1 = sel0.column(1);
181    let ymin = col1.min().unwrap();
182    let ymax = col1.max().unwrap();
183    let xdel = xmax - xmin;
184    let ydel = ymax - ymin;
185    let lbar = (xdel + ydel) / 2.0;
186    let ftol = ftol.unwrap_or(&5.0e-14);
187    let feps = ftol * lbar;
188    // Mask to filter vertices
189    let mask = vert
190        .outer_iter()
191        .map(|row| {
192            row[0] >= xmin - feps
193                && row[1] >= ymin - feps
194                && row[0] <= xmax + feps
195                && row[1] <= ymax + feps
196        })
197        .collect::<Vec<bool>>();
198    let selected_indices: Vec<usize> = mask
199        .iter()
200        .enumerate()
201        .filter_map(|(i, &m)| if m { Some(i) } else { None })
202        .collect();
203    let mut vert: Array2<f64> = vert.select(Axis(0), &selected_indices);
204    if vert.is_empty() {
205        return (STAT, BNDS);
206    }
207
208    // repeated code? xdel, ydel, lbar above were already computed from the subset
209    // let xdel = vert.slice(s![.., 0]).max().unwrap() - vert.slice(s![.., 0]).min().unwrap();
210    // let ydel = vert.slice(s![.., 1]).max().unwrap() - vert.slice(s![.., 1]).min().unwrap();
211    // let lbar = (xdel + ydel) / 2.0;
212
213    // flip to ensure y-axis is the `long` axis
214    let mut node = node.clone();
215
216    if xdel > ydel {
217        vert.invert_axis(Axis(1));
218        node.invert_axis(Axis(1));
219    };
220
221    // Sort edges by y-value
222    let swap = edge
223        .column(1)
224        .iter()
225        .zip(edge.column(0))
226        .map(|(&y1, &y0)| node[(y1, 1)] < node[(y0, 1)])
227        .collect::<Vec<bool>>();
228
229    let mut edge = edge.clone();
230
231    for (i, &swap) in swap.iter().enumerate() {
232        if swap {
233            edge.swap([i, 0], [i, 1]);
234        }
235    }
236    let (stat, bnds) = _inpoly(&vert, &node, &edge, ftol, &lbar);
237    for (i, &si) in selected_indices.iter().enumerate() {
238        STAT[si] = stat[i];
239        BNDS[si] = bnds[i];
240    }
241    (STAT, BNDS)
242}
243
244//
245#[cfg(test)]
246mod tests {
247
248    use crate::inpoly2;
249    use gnuplot::Caption;
250    use gnuplot::Color;
251    use gnuplot::Figure;
252    use gnuplot::PointSymbol;
253    use meshgridrs;
254    use ndarray::array;
255    use ndarray::stack;
256    use ndarray::Array;
257    use ndarray::Axis;
258
259    #[test]
260    fn example_1() {
261        // node = np.array([
262        //         [4, 0], [8, 4], [4, 8], [0, 4], [3, 3],
263        //         [5, 3], [5, 5], [3, 5]])
264        let node = array![
265            [4.0, 0.0],
266            [8.0, 4.0],
267            [4.0, 8.0],
268            [0.0, 4.0],
269            [3.0, 3.0],
270            [5.0, 3.0],
271            [5.0, 5.0],
272            [3.0, 5.0]
273        ];
274        // edge = np.array([
275        //     [0, 1], [1, 2], [2, 3], [3, 0], [4, 5],
276        //     [5, 6], [6, 7], [7, 4]]);
277        let edge = array![
278            [0, 1],
279            [1, 2],
280            [2, 3],
281            [3, 0],
282            [4, 5],
283            [5, 6],
284            [6, 7],
285            [7, 4]
286        ];
287        // xpos, ypos = np.meshgrid(
288        //     np.linspace(-1, 9, 51), np.linspace(-1, 9, 51))
289        let x = Array::linspace(-1., 9., 51);
290        let y = Array::linspace(-1., 9., 51);
291        let xi = vec![x, y];
292        let grids = meshgridrs::meshgrid(&xi, meshgridrs::Indexing::Xy).unwrap();
293        let xpos = &grids[0];
294        let ypos = &grids[1];
295        let xpos_flat = xpos.to_shape((xpos.len(), 1)).unwrap();
296        let ypos_flat = ypos.to_shape((ypos.len(), 1)).unwrap();
297        let points = stack![Axis(1), xpos_flat, ypos_flat].remove_axis(Axis(2));
298        let (inpoly, onedge) = inpoly2(&points, &node, Some(&edge), None);
299        let points_in: Vec<(f64, f64)> = points
300            .outer_iter()
301            .zip(&inpoly)
302            .filter(|&(_, &in_set)| in_set)
303            .map(|(row, _)| (row[0], row[1]))
304            .collect();
305
306        let points_out: Vec<(f64, f64)> = points
307            .outer_iter()
308            .zip(&inpoly)
309            .filter(|&(_, &in_set)| !in_set)
310            .map(|(row, _)| (row[0], row[1]))
311            .collect();
312
313        let points_on: Vec<(f64, f64)> = points
314            .outer_iter()
315            .zip(&onedge)
316            .filter(|&(_, &on_set)| on_set)
317            .map(|(row, _)| (row[0], row[1]))
318            .collect();
319
320        let x_in: Vec<f64> = points_in.iter().map(|&(x, _)| x).collect();
321        let y_in: Vec<f64> = points_in.iter().map(|&(_, y)| y).collect();
322
323        let x_out: Vec<f64> = points_out.iter().map(|&(x, _)| x).collect();
324        let y_out: Vec<f64> = points_out.iter().map(|&(_, y)| y).collect();
325
326        let x_on: Vec<f64> = points_on.iter().map(|&(x, _)| x).collect();
327        let y_on: Vec<f64> = points_on.iter().map(|&(_, y)| y).collect();
328
329        let mut fg = Figure::new();
330
331        {
332            let axes = fg.axes2d();
333            axes
334                //.set_aspect_ratio(1.0)
335                .points(
336                    &x_in,
337                    &y_in,
338                    &[Caption("IN==1"), Color("blue"), PointSymbol('O')],
339                )
340                .points(
341                    &x_out,
342                    &y_out,
343                    &[Caption("IN==0"), Color("red"), PointSymbol('O')],
344                )
345                .points(
346                    &x_on,
347                    &y_on,
348                    &[Caption("ON==1"), Color("magenta"), PointSymbol('S')],
349                );
350        }
351        fg.show().unwrap();
352    }
353}
inpoly/lib.rs

inpoly/
lib.rs
