FZJ-INM1-BDA / siibra-python

Software interfaces for interacting with brain atlases - Python client
Apache License 2.0
48 stars 10 forks source link

Regions extracted using same bounding box from BigBrain histology and isocortex segmentation are not aligned #11

Closed marcenko closed 3 years ago

marcenko commented 3 years ago

Given a bounding box, I tried to extract a region from the BigBrain histological space and the corresponding isocortex segmentation. Since I used the same bounding box for both, I would expect the resulting volumes to be aligned (up to the precision of the cortex segmentation), but they are not. This is the code I used:

import numpy as np
import brainscapes as bs
from nilearn import plotting
from cloudvolume import Bbox

bb_img = bs.bigbrain.BigBrainVolume(bs.spaces.BIG_BRAIN_HISTOLOGY.url)

# Define two points in BigBrain space
p0 = np.array((-3.979, -61.256, 3.906))
p1 = np.array((5.863, -55.356, -2.487))

def mm_to_vox(p0, p1, img, mip):
    # Bounding box needs to be defined in voxel space, so we need to apply the inverse affine matrix to the points
    inv_aff = np.linalg.inv(img.affine(mip))
    p0_vox = np.dot(inv_aff[:3, :3], p0) + inv_aff[:3, -1]
    p1_vox = np.dot(inv_aff[:3, :3], p1) + inv_aff[:3, -1]
    return p0_vox, p1_vox

# Read region from BigBrain
mip = 0
p0_vox, p1_vox = mm_to_vox(p0, p1, img=bb_img, mip=mip)
# Define bounding box
bbox = Bbox(p0_vox, p1_vox)
img = bb_img.Image(clip=bbox, mip=mip, force=True)

mask_url = "https://neuroglancer.humanbrainproject.eu/precomputed/BigBrainRelease.2015/classif/"
mask_img = bs.bigbrain.BigBrainVolume(mask_url)

# Bounding box needs to be redefined, as mask_img has different voxel space
p0_vox, p1_vox = mm_to_vox(p0, p1, img=mask_img, mip=mip)
# Define bounding box
bbox = Bbox(p0_vox, p1_vox)
mask = mask_img.Image(clip=bbox, mip=mip, force=True)

# Plot on top of each other. I would expect these images to be aligned, since I used the same bounding box, but they are not.
plotting.view_img(mask,
                  bg_img=img,
                  cmap="gray",
                  vim=0,
                  vmax=255,
                  resampling_interpolation="nearest",
                  symmetric_cmap=False)

I could imagine that there is an issue with how the bounding box coordiantes are rounded internally to match to the voxel space.

dickscheid commented 3 years ago

After reviewing this issue in detail with @ch-schiffer , we verified that this is an error in the metadata of the custom neuroglancer URL given in the bug description. The error will not occur when using curated siibra templates and masks. we will cleanup unused metadata on the neuroglancer sites to avoid this pitfall in teh future.