intersection_candidates_with_other_tree() does not find all candidates

[JosiahParry]

I was trying to give geo-index a spin with the use case of finding self-intersecting geometries (contiguity). At minimum, we expect each geometry to be an intersection candidate with itself. Using geo-index only 1 intersection candidate is found. When compared to rstar, all candidates are found.

Repro below. Please excuse any horrendous use of geoarrow 🙃

use geo::BoundingRect;
use geo::Polygon;
use geo_index::rtree::RTreeIndex;
use geoarrow::algorithm::{geo_index::RTree, native::Concatenate};
use geoarrow::array::{AsChunkedGeometryArray, PolygonArray};
use geoarrow::io::geojson::read_geojson;
use geoarrow::{trait_::GeometryArrayAccessor, GeometryArrayTrait};
use rstar::primitives::GeomWithData;
use rstar::{primitives::Rectangle, AABB};

// Find tree self-intersection canddiates using rstar
fn geo_contiguity(geom: Vec<Polygon>) {
    let to_insert = geom
        .iter()
        .enumerate()
        .map(|(i, gi)| {
            let rect = gi.bounding_rect().unwrap();
            let aabb =
                AABB::from_corners([rect.min().x, rect.min().y], [rect.max().x, rect.max().y]);

            GeomWithData::new(Rectangle::from_aabb(aabb), i)
        })
        .collect::<Vec<_>>();

    let tree = rstar::RTree::bulk_load(to_insert);
    let candidates = tree.intersection_candidates_with_other_tree(&tree);

    println!("rstar candidates:");
    for (left_idx, right_idx) in candidates {
        println!("{:?} {:?}", left_idx.data, right_idx.data);
    }
}

// Find tree self-intersection canddiates using geo-index
fn contiguity(geoms: &PolygonArray<i32>, node_size: usize) -> Vec<(usize, usize, f64)> {
    let tree = geoms.create_rtree_with_node_size(node_size);
    let mut _neighbors = Vec::new();

    let candidates = tree.intersection_candidates_with_other_tree(&tree);

    println!("geo-index candidates:");
    for (left_idx, right_idx) in candidates {
        println!("{:?} {:?}", left_idx, right_idx);
    }

    _neighbors
}

fn main() {
    let file = std::fs::File::open("src/guerry.geojson").unwrap();
    let reader = std::io::BufReader::new(file);

    // read from file
    let geoms = read_geojson(reader, None).unwrap();
    let polys = geoms
        .geometry()
        .unwrap()
        .as_ref()
        .as_polygon()
        .concatenate()
        .unwrap();

    println!("There are {} polygons", polys.len());
    contiguity(&polys, 10);

    let geo_polys = polys.iter_geo_values().collect::<Vec<_>>();
    geo_contiguity(geo_polys);
}

Output:

There are 116 polygons
geo-index candidates:
520 520
rstar candidates:
39 39
39 85
39 50
39 84
85 39
85 85
85 50
85 84
50 39
50 85
50 50
50 84
84 39
84 85
84 50
84 84
39 38
39 86
39 87
85 38
85 86
85 87
84 86
84 87
50 41
39 102
39 57
50 27
50 57
38 39
38 85
86 39
86 85
86 84
87 39
87 85
87 84
7 7
7 38
38 7
38 38
86 86
86 87
87 86
87 87
7 9
38 9
38 101
38 102
9 7
9 38
101 38
9 9
9 101
101 9
101 101
101 102
41 50
14 14
14 16
14 41
40 40
40 41
16 14
16 16
16 41
41 14
41 40
41 16
41 41
14 113
14 27
16 27
41 27
41 57
16 18
16 17
16 111
14 112
14 99
16 112
16 99
16 15
102 39
27 50
57 39
57 50
102 38
102 101
113 14
27 14
27 16
27 41
57 41
102 102
102 57
113 113
113 27
27 113
27 27
27 57
57 102
57 27
57 57
102 56
113 20
27 20
27 56
57 56
113 112
113 26
113 45
27 26
20 113
20 27
56 102
56 27
56 57
20 20
20 56
56 20
56 56
20 26
18 16
17 16
111 16
18 18
17 17
111 111
111 99
112 14
112 16
99 14
99 16
15 16
112 113
99 111
112 112
112 99
99 112
99 99
15 15
112 45
112 46
99 46
26 113
26 27
45 113
26 20
45 112
46 112
46 99
26 26
26 45
45 26
45 45
45 46
46 45
46 46
7 88
38 10
38 88
9 42
9 88
101 42
101 10
102 10
20 10
56 10
20 13
56 13
20 83
26 83
26 2
26 19
45 19
111 110
111 54
111 59
99 54
112 59
99 59
46 94
46 59
112 51
45 51
46 51
10 38
88 7
88 38
42 9
42 101
10 101
88 9
10 102
10 20
10 56
13 20
13 56
83 20
83 26
2 26
19 26
19 45
37 37
37 42
37 10
42 37
42 42
42 10
10 37
10 42
10 10
88 88
37 108
37 11
37 58
10 53
10 58
10 13
37 29
37 5
108 37
11 37
105 105
105 106
106 105
106 106
106 108
106 11
108 106
108 108
108 11
11 106
11 108
11 11
106 103
106 104
108 107
108 29
108 5
106 3
108 3
108 4
11 3
103 106
104 106
103 103
104 104
53 10
58 37
58 10
13 10
52 52
52 53
53 52
53 53
53 58
53 13
58 53
58 58
58 13
13 53
13 58
13 13
52 47
52 29
52 5
53 5
58 5
52 83
52 2
53 83
13 83
52 91
52 93
52 0
29 37
5 37
107 108
29 108
5 108
47 52
29 52
5 52
5 53
5 58
107 107
107 29
47 47
47 29
47 5
29 107
29 47
29 29
29 5
5 47
5 29
5 5
47 3
47 4
29 3
29 4
5 83
47 91
47 0
3 106
3 108
3 11
4 108
3 47
3 29
4 47
4 29
3 3
3 4
4 3
4 4
83 52
83 53
83 13
2 52
83 5
83 83
83 2
2 83
2 2
2 19
19 2
19 19
83 93
2 93
2 76
19 76
91 52
93 52
91 47
93 83
93 2
76 2
76 19
91 91
91 93
93 91
93 93
93 76
76 93
76 76
91 0
93 49
93 0
0 52
0 47
49 93
0 91
0 93
28 28
28 49
49 28
49 49
49 0
0 49
0 0
19 115
19 55
19 51
93 22
76 22
76 115
76 55
93 21
76 21
93 48
49 22
28 92
28 63
49 92
28 90
49 48
110 111
54 111
54 99
59 111
59 112
59 99
94 46
59 46
69 69
70 70
74 74
73 73
73 54
73 71
73 25
73 36
73 43
73 72
73 68
54 73
71 73
110 110
110 54
109 109
54 110
54 54
71 71
54 59
54 43
54 64
59 54
94 94
94 59
59 94
59 59
59 43
94 81
94 64
59 64
25 73
36 73
25 25
25 36
36 25
36 36
33 33
34 34
25 43
25 24
25 44
25 23
43 73
72 73
68 73
43 54
43 59
43 25
61 61
61 43
43 61
43 43
72 72
68 68
61 81
61 64
43 64
61 60
43 60
43 44
61 12
64 54
81 94
64 94
64 59
81 61
64 61
64 43
81 81
81 64
64 81
64 64
81 30
81 12
24 25
24 24
32 32
35 35
44 25
23 25
60 61
60 43
44 43
60 60
44 44
23 23
12 61
30 81
12 81
30 30
30 96
30 12
96 30
96 96
96 12
12 30
12 96
12 12
94 51
94 31
81 31
30 98
30 31
30 80
96 80
96 82
96 100
51 112
51 45
51 46
115 19
55 19
51 19
22 93
22 76
115 76
55 76
22 49
21 93
21 76
92 28
92 49
63 28
48 93
90 28
48 49
51 94
31 94
31 81
98 30
31 30
80 30
80 96
82 96
100 96
22 22
22 115
115 22
115 115
115 55
55 115
55 55
55 51
51 55
51 51
22 92
22 114
22 63
22 21
115 97
55 97
55 98
55 31
51 31
22 8
115 8
92 22
114 22
63 22
21 22
92 92
92 114
92 63
114 92
114 114
114 63
63 92
63 114
63 63
21 21
114 90
114 66
114 67
63 66
63 67
63 62
63 8
114 89
63 75
90 114
90 90
48 48
90 66
90 89
97 115
97 55
98 55
31 55
31 51
97 97
97 95
97 98
95 97
95 95
95 98
98 97
98 95
98 98
98 31
31 98
31 31
97 62
97 8
97 80
98 80
97 1
8 22
8 115
66 114
66 63
67 114
67 63
62 63
8 63
66 90
62 97
8 97
66 66
66 67
67 66
67 67
67 62
67 8
62 67
62 62
62 8
8 67
8 62
8 8
66 89
66 75
62 1
66 65
67 65
67 6
62 6
89 114
75 63
89 90
89 66
75 66
89 89
89 75
75 89
75 75
80 97
80 98
80 80
80 100
82 82
82 100
100 80
100 82
100 100
80 79
80 1
82 1
82 77
100 1
100 77
82 78
100 78
1 97
1 62
79 80
1 80
1 82
1 100
77 82
77 100
79 79
79 1
1 79
1 1
1 77
77 1
77 77
1 6
77 6
77 78
65 66
65 67
6 67
6 62
78 82
78 100
6 1
6 77
78 77
65 65
6 6
78 78

[JosiahParry]

I forgot to add the geojson!
guerry.geojson.zip

[kylebarron]

Thanks! I'll take a look when I have time. The search method is directly ported from upstream, and so that should be reliable (it has matched shapely/GEOS). The intersection candidates is my own addition.

Contributions welcome 😉 . I can describe/document the tree structure better if you're interested! (I actually have a half-written blog post on flatbush that I never finished)

[JosiahParry]

I can confirm that the search() method works well! I'm able to identify polygon queen contiguity with a function like this.

pub fn search_tree(left: OwnedRTree<f64>, targets: PolygonArray<i32>) -> Vec<(usize, usize, f64)> {
    let mut neighbors = Vec::new();
    let _ = targets.iter().enumerate().for_each(|(i, target)| {
        let geom = target.unwrap();
        let (min, max) = bounding_rect_polygon(&geom);
        let mut candidates = left.search(min[0], min[1], max[0], max[1]);
        candidates.sort();
        for c in candidates {
            if geom.to_geo().intersects(&targets.get(c).unwrap().to_geo()) {
                neighbors.push((i, c, 1.0));
            }
        }
    });
    neighbors
}

I tried toying around with the tree traversal but im a bit out of my depth when it comes to recursion and tree searching. I'll try and find some more time to maybe pluck around in the source code at a later point

pub fn two_tree_candidates(left: OwnedRTree<f64>, right: OwnedRTree<f64>) -> Vec<Vec<usize>> {
    let res = right
        .as_ref()
        .root()
        .children()
        .flat_map(|node| recursive_search(&node, &left))
        .collect::<Vec<_>>();
    res
}

fn recursive_search<'a>(
    node: &'a Node<'a, f64, RTreeRef<'a, f64>>,
    left: &OwnedRTree<f64>,
) -> Vec<Vec<usize>> {
    if node.is_leaf() {
        let candidate = left.search(node.min_x(), node.min_y(), node.max_x(), node.max_y());
        vec![candidate]
    } else if node.is_parent() {
        node.children()
            .flat_map(|child| recursive_search(&child, left))
            .collect()
    } else {
        Vec::new()
    }
}

[kylebarron]

I can confirm that the search() method works well!

Indeed, there are never any bugs in @ mourner's code 😉

I think the gist of my tree traversal is sound. It's probably just missing a check somewhere.

[kylebarron]

I put together a small derived repro in #51. Run with cargo run -- --debug. It's simplest to remove geoarrow from the picture here.

One of the issues is that

geo-index/src/rtree/traversal.rs

Line 56 in 159f2db

self.index >= self.tree.num_items() * 4

is backwards and should be < instead. So that meant that the iterator was never recursing:

geo-index/src/rtree/traversal.rs

Line 190 in 159f2db

(true, true) => return Some((left.index, right.index)),

Fixing that gives endless recursion, so I still have another issue somewhere. It's probably that I'm mixing parent and child intersections somewhere