searching for neighbours in the nth-ring around the center cell

[keewis]

Layer.neighbours currently allows searching for the immediate neighbours of a cell. However, I need to be able to query the n-th ring of neighbours around that cell (the aim is to create a space-domain convolution algorithm based on sparse matrices).

As far as I can tell, this should be possible by going integer offsets from the (face, x, y) (I think the code in cdshealpix calls this d0h, i, and j?) representation of the input cell:

let (new_face, new_x, new_y) = account_for_face_boundary(f, x + x_offset, y + y_offset);

assuming that account_for_face_boundary is a function that recalculates f, x, and y if any of x or y is outside the range of [0, nside[... this is similar to what happens in edge_cell_neighbours.

So I was wondering whether you'd be open to adding a function that allows searching for neighbours within a k-ring around the center cell (somewhat like h3ronpy's grid_disk)?

(As a warning, I'm working in the geo-sciences, and will probably open a couple more requests like this. If you think any of these are out-of-scope for this project I'd be fine with moving those to a separate, geo-specific crate, but would then need access to lower-level functions like decode or bits_2_hash)

[fxpineau]

@keewis, thank you for having open an issue.
This message to let you know that I am currently implementing the feature you are describing.

[fxpineau]

@keewis, done, but be aware of "strange" behaviors near base cell corners having only 2 base cell neighbours (instead of 3).
You may also have a look at the ring_coverage_approx method (I have to fix the doc: we see the copy/paste from cone_coverage_approx...)

Yes:

• d0h stands for depth 0 hash, i.e. the base cell number (in [0, 12[).
• i is the pixel index inside the base cell, in [0, nide[, along the South base cell vertex to East base cell vertex axis
• j is the pixel index inside the base cell, in [0, nide[, along the South base cell vertex to West base cell vertex axis

[fxpineau]

A few images to illustrate normal + corner cases:

[keewis]

Thanks, that was much quicker than I had hoped! What's the difference between ring_coverage_approx and ring_coverage_approx_custom?

A few images to illustrate normal + corner cases:

Yes, my own implementation (healpy-based, so didn't want to point you to that before you wrote it because of the licensing) also has that issue. For my application (the space-domain convolution), this tends to not matter too much, but it is definitely something to look out for.

[fxpineau]

Having a look at:

n_neighbours = (2 * ring + 1) ** 2

I understand that you expect all cells inside a square area of half size k cells and not only cells in the ring
of "distance" k (the internal border of the square area), it it right?
So far I implemented the latter, but it is easy to implement the first option (I have to replace 4 for loops by 2 nested for loops). Do you need both?

The ring_coverage_approx_custom method returns all cells in a given annulus defined by 2 radii (the internal radius and the external radius, both in angular distances).
If you want all cells inside a given angular radius, you may call the cone_coverage_approx_custom methods.
For a square shape on the sky, you may have a look at the box_coverage method (with a=b).
For visual examples, see the doc of cds-healpix-python, a python wrapper for the Rust library.

You may also have a look at the MOC libraries: the rust one and the MOCPy python wrapper, see a visual example in the MOCPy doc here.

The current test to know if a cell intersects a Cone (or a Ring) is based on a quick approximation leading to false positives.
With custom, you perform the computation at a lower depth before degrading the result at the requested depth.
It lowers the number of false positives at the cost of a larger computation time.

I still have to port in Rust the exact cell/cone intersection code I developed in Java.

[keewis]

So far I implemented the latter, but it is easy to implement the first option (I have to replace 4 for loops by 2 nested for loops). Do you need both?

No, I only need the one that returns every cell contained by the (discrete) ring around the center (so maybe it should have been called cone instead? Not sure, here my partial knowledge of healpix shows). In my own implementation this also returns the original cell id as the 0th ring, but if you don't want to also do that I can easily prepend that using the iterators.

[fxpineau]

Ok. I may call such a method k_neigbours or k_kneighbour_matrix, ... ?
Do you need the neighbours cells values to be always returned in the same order in a (2k+1)^2 matrix of Option<u64>
or a simple list of u64 (in any order) is ok?
(The latter will not change my code structure while the first option requires to re-work the code, testing for each pixel)

[keewis]

I've been returning ring by ring from the innermost (0th) to the outermost (kth) ring as a flattened list (the order within the ring is not important), but for my particular use case I only require the center cell to be first.

What do you return if there's only 7 neighbours for a corner cell? -1 like healpy? (won't work because this is u64, this would have to be the maximum value of u64) Or was that what you were referring to with Option<u64>?

[fxpineau]

The neighbours method returns a MainWindMap which has several methods such as get(MainWind) or entries (I may also have implemented IntoIter to get an Iterator). The get method returns an Option.

If you do not need the specific value for a given (x, y) position in the 2D matrix, then no need to use Options.

[keewis]

It would probably be simpler to still (have the option to) return Vec<Option<u64>> since that makes integration with cds-healpix-python easier (numpy arrays can't be ragged so the missing values will have to be filled, and either way it's nice to have some ordering in the flattened list).