From ea6d37dc239c16ff0d97d5fbc14a97c7443c6e46 Mon Sep 17 00:00:00 2001
From: Ross Markello <rossmarkello@gmail.com>
Date: Tue, 8 Jun 2021 14:16:31 -0400
Subject: [PATCH 1/2] [REF] Handle atlases w/decimal vox sizes

Adds improved handling of atlases with non-integer voxel sizes
---
 abagen/matching.py | 13 ++++++++-----
 1 file changed, 8 insertions(+), 5 deletions(-)

diff --git a/abagen/matching.py b/abagen/matching.py
index 0992a94..7666ea0 100644
--- a/abagen/matching.py
+++ b/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'

From 815257e9faa3a269ff1c617b788938aacaa4769c Mon Sep 17 00:00:00 2001
From: Ross Markello <rossmarkello@gmail.com>
Date: Thu, 10 Jun 2021 09:14:33 -0400
Subject: [PATCH 2/2] [TEST] Test for nonint voxels

---
 abagen/tests/test_matching.py | 19 +++++++++++++++++--
 1 file changed, 17 insertions(+), 2 deletions(-)

diff --git a/abagen/tests/test_matching.py b/abagen/tests/test_matching.py
index 64ab194..f6927ec 100644
--- a/abagen/tests/test_matching.py
+++ b/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