TEST: Centaurz2

.gitignore

@@ -1,6 +1,7 @@
 *.pyc
 MANIFEST
 pyxnat.egg-info/
+bbrc_pyxnat.egg-info/
 build/
 dist/
 *.DS_Store*

pyxnat/core/derivatives/3dasl.py

@@ -0,0 +1,19 @@
+from .basil import volumes as vols
+
+XNAT_RESOURCE_NAMES = ['3DASL_QUANTIFICATION']
+
+
+def volumes(self):
+    """Partial volume segmentations for CSF, GM, WM (fsl_anat)"""
+    return vols(self)
+
+
+def perfusion(self):
+    """Perfusion mean values and summary statistics"""
+    import pandas as pd
+    from io import StringIO
+
+    f = self.file('quantification_results.csv')
+    content = self._intf.get(f.attributes()['URI']).content.decode('utf-8')
+    df = pd.read_csv(StringIO(content))
+    return df

pyxnat/core/derivatives/__init__.py

@@ -8,25 +8,26 @@
 # aseg = resource.aseg()
 # and thus get access to FreeSurfer measurements.
 #
-# Another example, with the ASHS pipeline for hippocampal subfield segmentation:
+# Another example, with the ASHS pipeline for hippocampal subfield
+# segmentation:
 #
 # resource = experiment.resource('ASHS')
 # volumes = resource.volumes()
 #
 # Again, this would only work provided that corresponding resources respect a
-# certain naming and structure on XNAT. Here in this present example, FreeSurfer
-# results are stored in resources called FREESURFER6 and the whole FreeSurfer
-# folder (named after the subject) is stored in the resource. Having access
-# to such additional functions would be conditioned by the existence of these
-# resources with proper matching structure.
+# certain naming and structure on XNAT. Here in this present example,
+# FreeSurfer results are stored in resources called FREESURFER6 and the whole
+# FreeSurfer folder (named after the subject) is stored in the resource.
+# Having access to such additional functions would be conditioned by the
+# existence of these resources with proper matching structure.
 #
 # Nevertheless, this mechanism has been implemented so as to get easily adapted
 # to local configurations, by editing/adding this very same folder.
 #
 # Adding a custom function can be done simply as follows.
 #
-# In this same folder (pyxnat/core/derivatives/), edit any existing file or add a new
-# one (filename does not matter):
+# In this same folder (pyxnat/core/derivatives/), edit any existing file or add
+# a new one (filename does not matter):
 #
 # Define XNAT_RESOURCE_NAME. This variable names the XNAT resource which needs
 # a custom function.
@@ -37,4 +38,3 @@
 # XNAT_RESOURCE_NAMES.
 #
 # An example is provided in pyxnat/core/derivatives/ashs.py
-# 

pyxnat/core/derivatives/alps.py

@@ -0,0 +1,34 @@
+XNAT_RESOURCE_NAME = 'ALPS'
+
+
+def alps(self):
+    """ Returns the average ALPS (Diffusivity Along the Perivascular Space) index
+    for both hemispheres, along with the left and right side values and supporting
+    metrics."""
+
+    import pandas as pd
+    from io import StringIO
+
+    f = self.file('alps.stat/alps.csv')
+    content = self._intf.get(f.attributes()['URI']).content.decode('utf-8')
+    df = pd.read_csv(StringIO(content))
+
+    # drop empty columns and sort table
+    df.drop(columns=['id', 'scanner'], inplace=True)
+    return df
+
+
+def fa_md_alps(self):
+    """Returns FA (Fractional Anisotropy) and MD (Mean Diffusivity) diffusion
+    metrics measured in the same ROIs used for the ALPS index, for both hemispheres."""
+
+    import pandas as pd
+    from io import StringIO
+
+    f = self.file('alps.stat/fa+md_alps.csv')
+    content = self._intf.get(f.attributes()['URI']).content.decode('utf-8')
+    df = pd.read_csv(StringIO(content))
+
+    # drop empty columns and sort table
+    df.drop(columns=['id', 'scanner'], inplace=True)
+    return df

pyxnat/core/derivatives/ashs.py

@@ -1,32 +1,33 @@
 XNAT_RESOURCE_NAME = 'ASHS'
 
+
 def volumes(self, mode='corr_nogray'):
     import pandas as pd
 
     f = list(self.files('*icv.txt'))[0]
     uri = f._uri
     res = self._intf.get(uri).text.split(' ')
 
-    f = list(self.files('*_left_%s_volumes.txt'%mode))[0]
+    f = list(self.files('*_left_%s_volumes.txt' % mode))[0]
     uri = f._uri
     resl = self._intf.get(uri).text.split('\n')
 
-    f = list(self.files('*_right_%s_volumes.txt'%mode))[0]
+    f = list(self.files('*_right_%s_volumes.txt' % mode))[0]
     uri = f._uri
     resr = self._intf.get(uri).text.split('\n')
 
     table = []
     for resx in [resl, resr]:
         for line in resx:
-                if line == '':
-                    continue
-                line = line.split(' ')
-                s = line[0]
-                side = line[1]
-                region = line[2]
-                i = int(line[-2])
-                m = float(line[-1])
-                table.append([s, side, region, i, m])
+            if line == '':
+                continue
+            line = line.split(' ')
+            s = line[0]
+            side = line[1]
+            region = line[2]
+            i = int(line[-2])
+            m = float(line[-1])
+            table.append([s, side, region, i, m])
 
     table.append([res[0], None, 'tiv', None, float(res[1].rstrip('\n'))])
 

pyxnat/core/derivatives/bamos.py

@@ -0,0 +1,224 @@
+import os
+import nibabel as nib
+import tempfile
+import numpy as np
+import pandas as pd
+
+XNAT_RESOURCE_NAME = 'BAMOS'
+
+BAMOS_LABELS = {0: 'background',
+                1: 'left frontal',
+                2: 'right frontal',
+                3: 'left parietal',
+                4: 'right parietal',
+                5: 'left occipital',
+                6: 'right occipital',
+                7: 'left temporal',
+                8: 'right temporal',
+                9: 'basal ganglia',
+                10: 'infratentorial'}
+
+
+def volume(self):
+    """ Returns the estimated volume of the identified lesions (i.e. voxels
+    with a value higher than 0.5."""
+
+    fd, fp = tempfile.mkstemp(suffix='.nii.gz')
+    os.close(fd)
+
+    f = list(self.files('CorrectLesion_*nii.gz'))[0]
+    f.get(fp)
+    d = nib.load(fp)
+    size = np.prod(d.header['pixdim'].tolist()[:4])
+    n = np.array(d.dataobj)
+    v = np.sum(n[n >= 0.5]) * size
+    os.remove(fp)
+    return v
+
+
+def n_lesions(self):
+    """ Returns the estimated number of lesions i.e. the # of found connected
+    components in Connect_WS3WT3WC1Lesion*_corr.nii.gz. """
+
+    fd, fp = tempfile.mkstemp(suffix='.nii.gz')
+    os.close(fd)
+
+    f = list(self.files('Connect_WS3WT3WC1Lesion_*_corr.nii.gz'))[0]
+    f.get(fp)
+    d = nib.load(fp)
+    n = len(np.unique(d.dataobj))
+
+    os.remove(fp)
+    return n
+
+
+def stats(self):
+    """ Collects descriptive statistics based on the segmented lesions of the
+    white matter, including volumes and number of lesions. A voxel is
+    considered as a part of a lesion if it has a value higher than 0.5."""
+    def _download_data_(self):
+        fd, fp = tempfile.mkstemp(suffix='.nii.gz')
+        os.close(fd)
+
+        f = list(self.files('Connect_WS3WT3WC1Lesion_*_corr.nii.gz'))[0]
+        f.get(fp)
+        cc = np.array(nib.load(fp).dataobj)
+        f = list(self.files('CorrectLesion_*nii.gz'))[0]
+        f.get(fp)
+        d = nib.load(fp)
+        size = np.prod(d.header['pixdim'].tolist()[:4])
+        les = np.array(d.dataobj)
+
+        f = list(self.files('Layers_*.nii.gz'))[0]
+        f.get(fp)
+        d1 = np.array(nib.load(fp).dataobj)
+
+        f = list(self.files('Lobes_*.nii.gz'))[0]
+        f.get(fp)
+        d2 = np.array(nib.load(fp).dataobj)
+        os.remove(fp)
+        return cc, d1, d2, les, size
+
+    def _roistats_from_map(m, atlas1, atlas2, func=np.mean):
+        import itertools
+        assert(m.shape == atlas1.shape)
+        assert(m.shape == atlas2.shape)
+        labels1 = list(np.unique(atlas1))
+        labels2 = list(np.unique(atlas2))
+
+        combinations = list(itertools.product(labels1, labels2))
+
+        mask = lambda i1, i2: (atlas1 == i1) & (atlas2 == i2)
+        label_values = dict([((i1, BAMOS_LABELS[i2]), func(m[mask(i1, i2)]))
+                             for (i1, i2) in combinations])
+        return label_values
+
+
+    cc, d1, d2, les, size = _download_data_(self)
+    d1[d1 == 5] = 4  # Merging layer 5 with layer 4
+
+    stats1 = _roistats_from_map(cc, d1, d2, func=lambda d: len(np.unique(d)))
+    stats2 = _roistats_from_map(les, d1, d2, func=lambda d: np.sum(d[d >= 0.5]))
+
+    df = [(d, r, val, stats2[(d, r)] * size) for (d, r), val in stats1.items()]
+    df = pd.DataFrame(df, columns=['depth', 'region', 'n', 'volume'])
+
+    return df
+
+
+def bullseye_plot(self, ax=None, figsize=(12, 8), segBold=[],
+                  measurement='volume', stats=None, cm=None):
+    """
+    Bullseye's representation of cerebral white matter hyperintensities
+    (adapted from https://matplotlib.org/stable/gallery/specialty_plots/leftventricle_bulleye.html)
+
+    Args:
+    `ax`          If None then a new Matplotlib figure is created
+    `figsize`     Figure size
+    `segBold`     Highlight some sections (if an integer is passed then will
+                   highlight the regions over the given nth percentile)
+    `measurement` Measurement to plot (default: `volume`)
+    `stats`       Can be used to plot precalculated stats. (default: None -
+                   stats will be calculated for the current resource)
+    `cm`          Color map (default: `matplotlib.pyplot.cm.YlOrRd`)
+
+    Reference:
+    C. Sudre et al., J Neuroradiol, 2018
+    """
+    from matplotlib import pyplot as plt
+    import matplotlib as mpl
+
+    num = 80   # sampling parameter
+    nreg = 10  # how many regions
+    ndep = 4   # how many layers
+    offset = 18 * np.pi/180  # in radians
+    ordered_labels = [3, 1, 2, 4, 8, 6, 10, 9, 5, 7]
+
+    def _ravel_stats_(stats, measurement='n'):
+        data = []
+        for d in range(1, ndep+1):
+            for label in ordered_labels:
+                q = 'depth == %s & region == "%s"' % (d, BAMOS_LABELS[label])
+                v = float(stats.query(q)[measurement])
+                data.append(v)
+        return data
+
+    if stats is None:
+        stats = self.stats()
+
+    layers = [int(e) for e in set(stats['depth'])][1:]  # remove 0 (background)
+
+    # Highlight regions > nth percentile
+    if isinstance(segBold, int):
+        pc = stats[measurement].quantile(q=segBold/100.0)
+        segBold = []
+        for _, r in stats.query('%s > @pc' % measurement).iterrows():
+            label = [k for (k, v) in BAMOS_LABELS.items() if v == r.region][0]
+            idx = ordered_labels.index(label)
+            segBold.append((int(r.depth)-1, idx))
+
+    data = _ravel_stats_(stats, measurement=measurement)
+
+    # Start plotting
+    data = np.array(data).ravel()
+    vlim = [data.min(), data.max()]
+    if cm is None:
+        cm = plt.cm.YlOrRd
+
+    axnone = ax is None
+    if axnone:
+        fig, ax = plt.subplots(subplot_kw=dict(projection='polar'),
+                               figsize=figsize)
+
+    theta = np.linspace(0, 2*np.pi, num*nreg)
+    r = np.linspace(0.2, 1, ndep+1)
+
+    # Draw circles
+    linewidth = 1
+    for i in range(r.shape[0]):
+        ax.plot(theta, np.repeat(r[i], theta.shape), '-k', lw=linewidth)
+
+    # Draw lines
+    for i in range(nreg):
+        theta_i = (i * 360.0/nreg) * np.pi/180 + offset
+        ax.plot([theta_i, theta_i], [r[0], 1], '-k', lw=linewidth)
+
+    # Paint areas
+    for depth in range(ndep):
+        r0 = r[depth:depth+2]
+        r0 = np.repeat(r0[:, np.newaxis], num, axis=1).T
+        for i in range(nreg):
+            theta0 = theta[i*num:(i+1)*num] + offset
+            theta0 = np.repeat(theta0[:, np.newaxis], 2, axis=1)
+            z = np.ones((num, 2)) * data[depth*nreg + i]
+            ax.pcolormesh(theta0, r0, z, vmin=vlim[0], vmax=vlim[1], cmap=cm)
+
+            # Highlight some regions
+            if (depth, i) in segBold:
+                ax.plot(theta0, r0, '-k', lw=linewidth+2)
+                ax.plot(theta0[0], [r[depth], r[depth+1]], '-k', lw=linewidth+1)
+                ax.plot(theta0[-1], [r[depth], r[depth+1]], '-k', lw=linewidth+1)
+
+    ax.set_ylim([0, 1])
+    ax.set_yticklabels(layers)
+
+    thetaticks = np.arange(0, 360, 36)
+    ordered_labels = [7, 3, 1, 2, 4, 8, 6, 10, 9, 5]
+
+    labels = [BAMOS_LABELS[e] for e in ordered_labels]
+    ax.set_thetagrids(thetaticks, labels=labels)
+    ax.tick_params(pad=12, labelsize=14, axis='both')
+    bbox = dict(boxstyle="round", ec="white", fc="white", alpha=0.5)
+    plt.setp(ax.get_xticklabels(), bbox=bbox)
+
+    if axnone:  # Add legend
+        cNorm = mpl.colors.Normalize(vmin=vlim[0], vmax=vlim[1])
+
+        ax = fig.add_axes([0.3, 0.04, 0.45, 0.05])
+        ticks = [vlim[0], 0, vlim[1]]
+        cb = mpl.colorbar.ColorbarBase(ax, cmap=cm, norm=cNorm,
+                                       orientation='horizontal', ticks=ticks)
+    plt.show()
+
+    if axnone:
+        return fig, ax

pyxnat/core/derivatives/bamos_arterial.py

@@ -0,0 +1,17 @@
+XNAT_RESOURCE_NAME = 'BAMOS_ARTERIAL'
+
+
+def stats(self):
+    """ Collects descriptive statistics based on the segmented lesions of the
+        white matter, including volumes and number of lesions per arterial
+        territory. A voxel is considered as a part of a lesion if it has a
+        value higher than 0.5."""
+    import pandas as pd
+    from io import StringIO
+
+    f = self.file('bamos_arterial_stats.csv')
+    uri = f._uri
+    res = self._intf.get(uri).text
+    text = StringIO(res)
+    df = pd.read_csv(text)
+    return df