Merge pull request #197 from rmarkello/fix/decimal_vox

[REF] Handle atlases w/decimal vox sizes

abagen/matching.py

@@ -5,7 +5,6 @@
 
 import warnings
 
-import nibabel as nib
 import numpy as np
 import pandas as pd
 from scipy.spatial import cKDTree, distance_matrix
@@ -55,8 +54,11 @@ def __init__(self, atlas, coords=None, *, triangles=None, atlas_info=None,
                               'constructor but `coords` is not None. Ignoring '
                               'supplied `coords` and using coordinates '
                               'derived from image.')
-            self._volumetric = nib.affines.voxel_sizes(affine)
             self._shape = atlas.shape
+            self._volumetric = tuple()
+            vox = affine[:-1, :-1][np.where(affine[:-1, :-1])]  # TODO: oblique
+            for vs, off, ndim in zip(vox, affine[:-1, -1], self._shape):
+                self._volumetric += (np.arange(off, off + (vs * ndim), vs),)
             nz = atlas.nonzero()
             atlas, coords = atlas[nz], transforms.ijk_to_xyz(np.c_[nz], affine)
         except TypeError:
@@ -246,9 +248,10 @@ def label_samples(self, annotation, tolerance=2):
 
         if self.volumetric:
             if self.group_atlas:
-                cols = ['mni_x', 'mni_y', 'mni_z']
-                vox_size = 1 / self._volumetric
-                samples[cols] = np.floor(samples[cols] * vox_size) / vox_size
+                # floor sample MNI coordinates to the grid of the atlas
+                for n, col in enumerate(['mni_x', 'mni_y', 'mni_z']):
+                    idx = np.sort(self._volumetric[n])
+                    samples[col] = idx[np.searchsorted(idx, samples[col]) - 1]
             labels = self._match_volume(samples, abs(tolerance))
         else:
             cortex = samples['structure'] == 'cortex'

abagen/tests/test_matching.py

@@ -8,7 +8,7 @@
 import pandas as pd
 import pytest
 
-from abagen import images, matching
+from abagen import images, matching, transforms
 
 
 def test_check_label():
@@ -63,7 +63,7 @@ def test_closest_centroids():
     assert np.allclose(out, 0)
 
 
-def test_AtlasTree(atlas, surface, testfiles):
+def test_AtlasTree(atlas, surface):
     # basic test data
     data = np.array([1, 1, 1, 2, 2, 2])
     coords = np.array([
@@ -166,3 +166,18 @@ def test_AtlasTree(atlas, surface, testfiles):
     with pytest.raises(ValueError):
         matching.AtlasTree(np.random.choice(10, size=(100,)),
                            coords=np.random.rand(99, 3))
+
+
+def test_nonint_voxels(atlas):
+    coord = [[-56.8, -50.6, 8.8]]
+    affine = np.zeros((4, 4))
+    affine[:-1, :-1] = np.eye(3) * 1.5
+    affine[:, -1] = nib.load(atlas['image']).affine[:, -1]
+    dataobj = np.zeros((100, 100, 100))
+    i, j, k = transforms.xyz_to_ijk(coord, affine).squeeze()
+    dataobj[i, j, k] = 1
+    newatl = nib.Nifti1Image(dataobj, affine)
+
+    tree = matching.AtlasTree(newatl, group_atlas=True)
+    assert tree.volumetric
+    assert tree.label_samples(coord, tolerance=0).loc[0, 'label'] == 1