Rounding scheme for neuroglancer fetching needs revision

[dickscheid]

When fetching an axis aligned flat volume of interest translating to exactly 1 voxel "thickness" in one dimension, the NeuroglancerScale fetching method rounds up to a 2 voxel array. In my example, fetching in bigbrain space with Bounding box from (36.58,-0.10,8.33)mm to (38.64,-0.08,11.71)mm in BigBrain microscopic template (histology) space retrieved an image array of shape (160, 2, 99).

This is critical, since several analysis workflows that I had implemented with siibra rely on the possibility to fetch a 2D patch from properly sized bounding boxes in physical coordinates, and would now confront me to resolve the ambiguity downstream.

The problem is caused by taking floorof the minpoint and ceil of the max point to compute the target shape as their difference, resulting in a doubling of integer rounding in the extreme case. The suggested change computes the rounded target shape as the ceil of the difference of the unrounded max- and midpoints instead, and adds it to the floored midpoint to obtain the rounded max point, therefore always only triggering a single rounding effect. Using the fix, the aforementioned bounding box yields an array of shape (160, 1, 98) after fetching, as expected.

The previous explicit additional clipping of voxel dimensions smaller than 1 seemed incorrect to me, and did not influence the result after the fix anyways, so I removed it.

The proposed change fixes the said problem, but it should be carefully thought through and tested on other examples. I suggest to define a range of standard and border cases of VOIs to double-check that fetching works as expected:

• perfect cube should result in perfect cube
• VOI dimensions equal or smaller to the volume's pixel size (resolution_mm) should result in flat dimensions (dim=1) of the flat array (that's the example reported above)
• VOIs larger than the volume dimensions should be clipped exactly to the max dimension of the volume
• ...

[codecov]

Codecov Report

Attention: Patch coverage is 50.00000% with 2 lines in your changes missing coverage. Please review.

Project coverage is 45.32%. Comparing base (f88a7d5) to head (dff053b).
Report is 97 commits behind head on main.

Files with missing lines	Patch %	Lines
siibra/volumes/providers/neuroglancer.py	50.00%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #642      +/-   ##
==========================================
- Coverage   45.34%   45.32%   -0.03%     
==========================================
  Files          75       75              
  Lines        7229     7228       -1     
==========================================
- Hits         3278     3276       -2     
- Misses       3951     3952       +1     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[xgui3783]

When fetching an axis aligned flat volume of interest translating to exactly 1 voxel "thickness" in one dimension, the NeuroglancerScale fetching method rounds up to a 2 voxel array. In my example, fetching in bigbrain space with Bounding box from (36.58,-0.10,8.33)mm to (38.64,-0.08,11.71)mm in BigBrain microscopic template (histology) space retrieved an image array of shape (160, 2, 99).

This is by design, no? We precisely do not want the user to deal in voxel space. So when bbox is provided, we try to satistify the user's request based on:

• user's request of resolution level
• the smallest collection of voxels that intersects with user's defined bbox

This will produce some weird edge cases - which is significantly magnified when user provides bbox with dimension similar or smaller than that of the voxel size.

In your case, if you definitely want a single voxel thickness, I would argue you would be fetching from (36.58,-0.09,8.33)mm to (38.64,-0.09,11.71)mm

I suggest to define a range of standard and border cases of VOIs to double-check that fetching works as expected:

The two guidelines above aimed to do just that (but perhaps poorly documented)

perfect cube should result in perfect cube

This is already like not possible, since the 1um slices cannot have perfect cube larger than 1 x 1 x 1

VOI dimensions equal or smaller to the volume's pixel size (resolution_mm) should result in flat dimensions (dim=1) of the flat array (that's the example reported above)

This will create weird edge cases where we now need to make a decision, if the bbox sits exactly in the middle of two voxels. We will have to at least create one more rule to disambiguate this edge case

---

edit:

I can't comment on the comment at the moment. I can take a look on monday.

[AhmetNSimsek]

(36.58,-0.10,8.33)mm to (38.64,-0.08,11.71)mm corrosponds to voxels bbox from (3330.127590551181, 3495.9708333333338, 5164.0913700787405) to (3170.4425511811023, 3494.9708333333338, 5066.7685354330715). While the difference for the y axis have exactly one voxel length, the question is defining how exactly siibra should handle data with 1 voxel length and intersecting with two voxels. Currently, siibra returns both of the voxels by taking the floor of the lower and ceil of the higher. One can use round for both to return exactly one voxel but this standard is not defined in siibra yet. This PR will be morphed to determine this standard.

[AhmetNSimsek]

removing dim increase is necessary. However, imo, the voxel selection calculation is correct. Please use debug mode to test the provided bounding boxes; if a bounding box touches more than one voxel, they are both selected. I could not reproduce a naturally occuring example that, say selects 2 voxels by floating errors.