Repository: AMYPAD/AmyPET (language: python)

Issue:
amypet.web requires `TKinter` module
Failed to run the web UI version of AmyPET. Looks like an additional package named `TKinter` is needed. `TKinter` is a Python interface to the Tcl/Tk toolkit. See trace below,

```
(amypet_dev) user@host:/.../tmp/amypet$ pip install amypet[web] --quiet

(amypet_dev) user@host:/.../tmp/amypet$ amypet.web
Traceback (most recent call last):
  File "/home/user/venvs/amypet_dev/bin/amypet.web", line 5, in <module>
    from amypet.web import main
  File "/home/user/venvs/amypet_dev/lib/python3.7/site-packages/amypet/web.py", line 3, in <module>
    from . import _backend_web
  File "/home/user/venvs/amypet_dev/lib/python3.7/site-packages/amypet/_backend_web.py", line 4, in <module>
    import tkinter as tk
ModuleNotFoundError: No module named 'tkinter'
```

I can confirm error can be suppressed by installing the missing module as follows,
```
apt-get install python3-tk
```

If such module is a requirement I would suggest noting it in the README file.
installation testing & docs
- [ ] windows (`conda`)
- [ ] linux (`conda`)
- [ ] linux (`apt`/`yum`)
- [ ] document in individual READMEs as well as one unified place, see DAG for:
  - https://github.com/AMYPAD/AmyPET
    - https://github.com/AMYPAD/NiftyPAD
    - https://github.com/NiftyPET/NiftyPET
      - https://github.com/NiftyPET/NIMPA
        - https://github.com/NiftyPET/NIPET
          - https://github.com/NiftyPET/NInst
        - https://github.com/AMYPAD/SPM12
        - https://github.com/AMYPAD/miutil
        - https://github.com/AMYPAD/CuVec
        - https://github.com/AMYPAD/NumCu
- [ ] p3: how does above DAG correspond to pipelines (i.e. what does centiloid require?)
- [ ] basic installation (default) and advanced installation (optional)

----

- required: Tk #23
- required: C++ compiler... Or make CuVec optional?
- optional (for which pipelines?): MATLAB #24
- optional (for which pipelines?): CUDA SDK


Discussion / hints:

Fix patch:
diff --git a/.github/workflows/comment-bot.yml b/.github/workflows/comment-bot.yml
index b44ee7ba..2a87fe42 100644
--- a/.github/workflows/comment-bot.yml
+++ b/.github/workflows/comment-bot.yml
@@ -1,25 +1,23 @@
 name: Comment Bot
 on:
-  issue_comment:
-    types: [created]
-  pull_request_review_comment:
-    types: [created]
+  issue_comment: {types: [created]}
+  pull_request_review_comment: {types: [created]}
 jobs:
   tag:  # /tag <tagname> <commit>
     if: startsWith(github.event.comment.body, '/tag ')
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
     - name: React Seen
-      uses: actions/github-script@v2
+      uses: actions/github-script@v6
       with:
         script: |
-          const perm = await github.repos.getCollaboratorPermissionLevel({
+          const perm = await github.rest.repos.getCollaboratorPermissionLevel({
             owner: context.repo.owner, repo: context.repo.repo,
             username: context.payload.comment.user.login})
           post = (context.eventName == "issue_comment"
-            ? github.reactions.createForIssueComment
-            : github.reactions.createForPullRequestReviewComment)
+            ? github.rest.reactions.createForIssueComment
+            : github.rest.reactions.createForPullRequestReviewComment)
           if (!["admin", "write"].includes(perm.data.permission)){
             post({
               owner: context.repo.owner, repo: context.repo.repo,
@@ -40,12 +38,12 @@ jobs:
         BODY: ${{ github.event.comment.body }}
         GITHUB_TOKEN: ${{ secrets.GH_TOKEN }}
     - name: React Success
-      uses: actions/github-script@v2
+      uses: actions/github-script@v6
       with:
         script: |
           post = (context.eventName == "issue_comment"
-            ? github.reactions.createForIssueComment
-            : github.reactions.createForPullRequestReviewComment)
+            ? github.rest.reactions.createForIssueComment
+            : github.rest.reactions.createForPullRequestReviewComment)
           post({
             owner: context.repo.owner, repo: context.repo.repo,
             comment_id: context.payload.comment.id, content: "rocket"})
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index b5f649e4..ab5c896f 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -2,7 +2,7 @@ default_language_version:
   python: python3
 repos:
 - repo: https://github.com/pre-commit/pre-commit-hooks
-  rev: v4.0.1
+  rev: v4.3.0
   hooks:
   - id: check-added-large-files
   - id: check-case-conflict
@@ -38,10 +38,11 @@ repos:
     - flake8-isort
     - flake8-string-format
 - repo: https://github.com/google/yapf
-  rev: v0.31.0
+  rev: v0.32.0
   hooks:
   - id: yapf
     args: [-i]
+    additional_dependencies: [toml]
 - repo: https://github.com/PyCQA/isort
   rev: 5.10.1
   hooks:
@@ -49,4 +50,4 @@ repos:
 - hooks:
   - id: mbeautify
   repo: https://github.com/AMYPAD/miutil
-  rev: v0.9.1
+  rev: v0.10.1
diff --git a/README.rst b/README.rst
index b133d7a1..e9e7c7a7 100644
--- a/README.rst
+++ b/README.rst
@@ -9,13 +9,11 @@ Amyloid imaging to prevent Alzheimer's Disease.
 Install
 -------
 
-Requires Python 3.6 or greater. Choose one of the following:
+Requires Python 3.6 or greater.
 
 .. code:: sh
 
-    pip install amypet[web]  # CLI and web user interface (Web UI) version (recommended)
-    pip install amypet       # CLI version
-    pip install amypet[gui]  # CLI and graphical user interface (GUI) version
+    pip install amypet
 
 For certain functionality (image trimming & resampling, DICOM to NIfTI converters) it may be required to also `pip install nimpa`.
 
@@ -27,7 +25,6 @@ Usage
 
     amypet.web     # Web UI version
     amypet --help  # CLI version
-    amypet.gui     # GUI version
 
 
 .. |Logo| image:: https://amypad.eu/wp-content/themes/AMYPAD/images/AMYPAD_Logo.jpg
diff --git a/amypet/__init__.py b/amypet/__init__.py
index f7b99085..ab91b3ee 100644
--- a/amypet/__init__.py
+++ b/amypet/__init__.py
@@ -20,3 +20,7 @@
         __licence__ = (Path(__file__).parent.parent / "LICENCE.md").read_text()
     except FileNotFoundError:
         __licence__ = "MPL-2.0"
+
+from .aux import *
+from .preproc import *
+from .suvr_tools import *
diff --git a/amypet/_backend_web.py b/amypet/_backend_web.py
index 7228c39e..e90cac14 100644
--- a/amypet/_backend_web.py
+++ b/amypet/_backend_web.py
@@ -1,7 +1,6 @@
 import logging
 import shlex
 import sys
-import tkinter as tk
 from argparse import (
     SUPPRESS,
     _HelpAction,
@@ -38,10 +37,6 @@
 
 {__licence__}"""}
 
-root = tk.Tk()
-root.withdraw()
-root.wm_attributes('-topmost', 1)
-
 
 class MyParser(BaseParser):
     def add_argument(self, *args, **kwargs):
@@ -88,18 +83,9 @@ def recurse(parser, key_prefix=""):
                             'name', NONE)
                     elif opt.widget == "DirChooser":
                         # https://github.com/streamlit/streamlit/issues/1019
-                        left, right = st.columns(2)
-                        with left:
-                            st.caption(opt.dest)
-                        with right:
-                            clicked = st.button("Browse directories", **kwargs)
-                        key = f'{key_prefix}{opt.dest}_val'
-                        if clicked:
-                            st.session_state[key] = tk.filedialog.askdirectory(
-                                master=root, initialdir=dflt) or dflt
-                        val = st.session_state.get(key, dflt)
-                        with left:
-                            st.write(f"`{val or '(blank)'}`")
+                        val = st.text_input(opt.dest, value=dflt, **kwargs)
+                        if val.startswith(prefix := "file://"):
+                            val = val[len(prefix):]
                     elif opt.widget == "IntegerField":
                         dflt = opt.default or 0
                         val = st.number_input(opt.dest, min_value=int(opt.widget_options['min']),
diff --git a/amypet/aux.py b/amypet/aux.py
new file mode 100644
index 00000000..6adb4ce8
--- /dev/null
+++ b/amypet/aux.py
@@ -0,0 +1,559 @@
+"""Auxiliary functions for the centiloid project"""
+__author__ = "Pawel Markiewicz"
+__copyright__ = "Copyright 2022"
+import pickle
+import re
+import urllib
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+from niftypet import nimpa
+from pkg_resources import resource_filename
+from scipy.stats import linregress
+
+# > regions used in CL project
+rvois = ['wc', 'cg', 'wcb', 'pns']
+
+# > folder name with the default conversion tables after calibration
+cl_anchor_fldr = Path(resource_filename(__name__, 'data/cl_anchor_tables'))
+cl_masks_fldr = Path(resource_filename(__name__, 'data/CL_masks'))
+amypet_dir = Path.home() / '.amypet'
+
+# > region full name strings for plots
+rvoi_str = {
+    'wc': 'WHOLE CEREBELLUM', 'cg': 'CEREBELLUM GM', 'wcb': 'WHOLE CEREBELLUM + BRAIN STEM',
+    'pns': 'PONS'}
+
+
+# ----------------------------------------------------------------------
+def get_atlas(atlas='aal', res=1, outpath=amypet_dir):
+    '''Get a brain atlas from `neuroparc` out of many available in MNI space.
+
+       Options:
+       - atlas:     one of many atlases, e.g., 'aal' which is the default;
+                    the number of atlases available can be extended by entering it
+                    in the dictionary. TODO: move the dictionary to `defs.py`.
+       - res:       the resolution of the atlas in mm.
+    '''
+    atlases = {
+        'ghpath': 'https://github.com/neurodata/neuroparc/raw/master/atlases/label/Human/',
+        'aal': 'AAL_space-MNI152NLin6_res-', 'aal_lbl': 'Anatomical-labels-csv/AAL.csv',
+        'schaefer200': 'Schaefer200_space-MNI152NLin6_res-',
+        'talairach': 'Talairach_space-MNI152NLin6_res-'}
+
+    # > get the GitHub data object
+    ghd = urllib.request.urlopen(atlases['ghpath'] + atlases[atlas] + f'{res}x{res}x{res}.nii.gz')
+    data = ghd.read()
+
+    # > save to file
+    nimpa.create_dir(outpath)
+
+    fatl = outpath / f'atlas-{atlas}_res-{res}mm.nii.gz'
+    with open(fatl, 'wb') as f:
+        f.write(data)
+
+    return fatl
+
+
+# ----------------------------------------------------------------------
+
+
+# ----------------------------------------------------------------------
+def im_check_pairs(fpets, fmris):
+    '''
+    checks visually image by image if the PET and MR have
+    the same orientation.
+    '''
+    fig, ax = plt.subplots(2, 2, figsize=(10, 10))
+    for fp, fm in zip(fpets, fmris):
+        print(f'{fp} : {fm}')
+        p = nimpa.getnii(fp)
+        m = nimpa.getnii(fm)
+
+        ps = p.shape
+        ms = m.shape
+
+        ax[0, 0].imshow(p[ps[0] // 2, ...], cmap='magma')
+        ax[1, 0].imshow(m[ms[0] // 2, ...], cmap='bone')
+
+        ax[0, 1].imshow(p[..., ps[2] // 2], cmap='magma')
+        ax[1, 1].imshow(m[..., ms[2] // 2], cmap='bone')
+
+        plt.draw()
+        plt.waitforbuttonpress(0)
+
+
+# ----------------------------------------------------------------------
+
+
+# ----------------------------------------------------------------------
+def identity_line(ax=None, ls='--', *args, **kwargs):
+    '''
+    plot identity line for any correlation analysis
+    needed for the CL calibration and beyond.
+    '''
+    # see: https://stackoverflow.com/q/22104256/3986320
+    ax = ax or plt.gca()
+    identity, = ax.plot([], [], ls=ls, *args, **kwargs)
+
+    def callback(axes):
+        low_x, high_x = ax.get_xlim()
+        low_y, high_y = ax.get_ylim()
+        low = min(low_x, low_y)
+        high = max(high_x, high_y)
+        identity.set_data([low, high], [low, high])
+
+    callback(ax)
+    ax.callbacks.connect('xlim_changed', callback)
+    ax.callbacks.connect('ylim_changed', callback)
+    return ax
+
+
+# ----------------------------------------------------------------------
+
+
+# ----------------------------------------------------------------------
+def get_clref(fxls):
+    '''
+    get the excel sheets with CL reference data from the CL project
+    excel sheet.
+    '''
+    import openpyxl
+
+    # > Define variable to load the workbook
+    wb = openpyxl.load_workbook(fxls)
+    ws = wb.active
+
+    pib_tbl = {'yc': {}, 'ad': {}}
+
+    # > cell offsets
+    ioff_ad = [5, 50]
+    ioff_yc = [51, 85]
+
+    sbj = [str(i.value) for i in ws['A']]
+    cg = [str(i.value) for i in ws['B']]
+    wc = [str(i.value) for i in ws['C']]
+    wcb = [str(i.value) for i in ws['D']]
+    pns = [str(i.value) for i in ws['E']]
+
+    ccg = [str(i.value) for i in ws['F']]
+    cwc = [str(i.value) for i in ws['G']]
+    cwcb = [str(i.value) for i in ws['H']]
+    cpns = [str(i.value) for i in ws['I']]
+
+    # > young control indices and SUVr/CL values
+    pib_tbl['yc']['id'] = np.array([int(i[3:]) for i in sbj[ioff_yc[0]:ioff_yc[1]]])
+
+    pib_tbl['yc']['suvr'] = {}
+    pib_tbl['yc']['suvr']['cg'] = np.array([float(i) for i in cg[ioff_yc[0]:ioff_yc[1]]])
+    pib_tbl['yc']['suvr']['wc'] = np.array([float(i) for i in wc[ioff_yc[0]:ioff_yc[1]]])
+    pib_tbl['yc']['suvr']['wcb'] = np.array([float(i) for i in wcb[ioff_yc[0]:ioff_yc[1]]])
+    pib_tbl['yc']['suvr']['pns'] = np.array([float(i) for i in pns[ioff_yc[0]:ioff_yc[1]]])
+
+    pib_tbl['yc']['cl'] = {}
+    pib_tbl['yc']['cl']['cg'] = np.array([float(i) for i in ccg[ioff_yc[0]:ioff_yc[1]]])
+    pib_tbl['yc']['cl']['wc'] = np.array([float(i) for i in cwc[ioff_yc[0]:ioff_yc[1]]])
+    pib_tbl['yc']['cl']['wcb'] = np.array([float(i) for i in cwcb[ioff_yc[0]:ioff_yc[1]]])
+    pib_tbl['yc']['cl']['pns'] = np.array([float(i) for i in cpns[ioff_yc[0]:ioff_yc[1]]])
+
+    # > AD indices and SUVr/CL
+    pib_tbl['ad']['id'] = np.array([int(i[3:]) for i in sbj[ioff_ad[0]:ioff_ad[1]]])
+    pib_tbl['ad']['suvr'] = {}
+    pib_tbl['ad']['suvr']['cg'] = np.array([float(i) for i in cg[ioff_ad[0]:ioff_ad[1]]])
+    pib_tbl['ad']['suvr']['wc'] = np.array([float(i) for i in wc[ioff_ad[0]:ioff_ad[1]]])
+    pib_tbl['ad']['suvr']['wcb'] = np.array([float(i) for i in wcb[ioff_ad[0]:ioff_ad[1]]])
+    pib_tbl['ad']['suvr']['pns'] = np.array([float(i) for i in pns[ioff_ad[0]:ioff_ad[1]]])
+
+    pib_tbl['ad']['cl'] = {}
+    pib_tbl['ad']['cl']['cg'] = np.array([float(i) for i in ccg[ioff_ad[0]:ioff_ad[1]]])
+    pib_tbl['ad']['cl']['wc'] = np.array([float(i) for i in cwc[ioff_ad[0]:ioff_ad[1]]])
+    pib_tbl['ad']['cl']['wcb'] = np.array([float(i) for i in cwcb[ioff_ad[0]:ioff_ad[1]]])
+    pib_tbl['ad']['cl']['pns'] = np.array([float(i) for i in cpns[ioff_ad[0]:ioff_ad[1]]])
+
+    return pib_tbl
+    # ----------------------------------------------------------------------
+
+
+def check_suvrs(suvr_yc, suvr_ad, refs):
+    ''' Check/QC the obtained SUVRs in `suvr_yc` and `suvr_ad `dictionaries
+        in relation to the reference in `refs` dictionary for the groups
+        ('ad' or 'yc').
+
+        Returns dictionary with the mean SUVrs and differences.
+    '''
+
+    # > prepare diff dictionary for the differences observed
+    # > in young controls (yc) and AD patients (ad)
+    diff = {'yc': {}, 'ad': {}}
+
+    # > initialise the means for each reference VOI
+    for rvoi in rvois:
+        diff['yc'][rvoi] = {'mean_ref': 0, 'mean': 0, 'N': 0}
+        diff['ad'][rvoi] = {'mean_ref': 0, 'mean': 0, 'N': 0}
+
+    def run_checks(suvr_dct, grp):
+
+        for rvoi in rvois:
+            print('========================================================')
+            for k in suvr_dct:
+
+                if grp == 'yc':
+                    idx = int(k[2:5])
+                elif grp == 'ad':
+                    idx = int(k[2:4])
+                else:
+                    raise ValueError('e> unknown group - only <yc>  or <ad> are accepted')
+
+                i = np.where(refs[grp]['id'] == idx)[0][0]
+
+                suvr = suvr_dct[k]['suvr'][rvoi]
+                suvr_ref = refs[grp]['suvr'][rvoi][i]
+                err = 100 * (suvr-suvr_ref) / suvr_ref
+
+                diff[grp][rvoi]['N'] += 1
+                diff[grp][rvoi]['mean_ref'] += suvr_ref
+                diff[grp][rvoi]['mean'] += suvr
+                diff[grp][rvoi][k] = {'suvr': suvr, 'ref': suvr_ref, 'err': err}
+
+                print(f'refvoi={rvoi}, indx> {idx}, suvr={suvr:.3f},'
+                      f' ref={suvr_ref:.3f}, error={err:.3f}%')
+
+            diff[grp][rvoi]['mean'] /= diff[grp][rvoi]['N']
+            diff[grp][rvoi]['mean_ref'] /= diff[grp][rvoi]['N']
+
+            # relative % mean difference
+            emean = diff[grp][rvoi]['mean']
+            rmean = diff[grp][rvoi]['mean_ref']
+            rmd = (emean-rmean) / rmean
+            diff[grp][rvoi]['mean_diff'] = rmd
+            print('--------------------------------------------------------')
+            print(f'> group % mean difference: {rmd:.3f}% (suvr={emean:.3f}, ref={rmean:.3f})')
+
+    print('========================================================')
+    print('YOUNG CONTROLS (YC)')
+    print('========================================================')
+    run_checks(suvr_yc, 'yc')
+
+    print('========================================================')
+    print('AD PATIENTS (AD)')
+    print('========================================================')
+    run_checks(suvr_ad, 'ad')
+
+    return diff
+
+
+# ----------------------------------------------------------------------
+
+
+# ----------------------------------------------------------------------
+def check_cls(suvr_yc, suvr_ad, diff, refs):
+    ''' check centiloid values relative to the reference in `refs`
+        dictionary, using SUVr dictionaries for YC and AD groups
+        `suvr_yc` and `suvr_ad`.  Also, update the `diff` dictionary
+        with the centiloid results and output it.
+    '''
+
+    for rvoi in rvois:
+        print('=====================================================')
+        suvr_0 = diff['yc'][rvoi]['mean']
+        suvr_100 = diff['ad'][rvoi]['mean']
+
+        grp = 'yc'
+        suvr_dct = suvr_yc
+
+        for k in suvr_dct:
+
+            idx = int(k[2:5])
+            i = np.where(refs[grp]['id'] == idx)[0][0]
+
+            cl = 100 * (suvr_dct[k]['suvr'][rvoi] - suvr_0) / (suvr_100-suvr_0)
+            cl_ref = refs[grp]['cl'][rvoi][i]
+
+            diff[grp][rvoi][k].update(cl=cl, cl_ref=cl_ref)
+
+            print(f'refvoi={rvoi}, indx> {idx}, cl={cl:.3f}, ref={cl_ref:.3f}')
+
+        grp = 'ad'
+        suvr_dct = suvr_ad
+
+        for k in suvr_dct:
+
+            idx = int(k[2:4])
+            i = np.where(refs[grp]['id'] == idx)[0][0]
+
+            cl = 100 * (suvr_dct[k]['suvr'][rvoi] - suvr_0) / (suvr_100-suvr_0)
+            cl_ref = refs[grp]['cl'][rvoi][i]
+
+            diff[grp][rvoi][k].update(cl=cl, cl_ref=cl_ref)
+
+            print(f'refvoi={rvoi}, indx> {idx}, cl={cl:.3f}, ref={cl_ref:.3f}')
+
+    def get_cls(diff, grp, rvoi, cl='cl'):
+        ''' Get centiloid values in a list for any group and reference VOI.
+        '''
+        return [diff[grp][rvoi][k][cl] for k in diff[grp][rvoi] if grp in k.lower()]
+
+    # ------------------------------------------------------------------
+    fig, ax = plt.subplots(2, 2, figsize=(10, 10))
+
+    for a, rvoi in enumerate(rvois):
+        i = a // 2
+        j = a % 2
+
+        # > combine CLs for both groups
+        cl = get_cls(diff, 'yc', rvoi, cl='cl')
+        cl += get_cls(diff, 'ad', rvoi, cl='cl')
+
+        # > combine reference CLs into one list
+        clref = get_cls(diff, 'yc', rvoi, cl='cl_ref')
+        clref += get_cls(diff, 'ad', rvoi, cl='cl_ref')
+
+        # > find the linear regression between the CL and reference CL
+        m, a, r, p, stderr = linregress(clref, cl)
+        r2 = r**2
+
+        # > plot the scatter plots for all 4 reference VOIs
+        ax[i, j].scatter(clref, cl, c='black')
+        identity_line(ax=ax[i, j], ls='--', c='b')
+        ax[i, j].text(60, 20, f'$y = {m:.4f}x + {a:.4f}$', fontsize=12)
+        ax[i, j].text(60, 10, f'$R^2={r2:.4f}$', fontsize=12)
+        ax[i, j].set_title(rvoi_str[rvoi])
+        ax[i, j].set_xlabel('CL ref')
+        ax[i, j].set_ylabel('CL AmyPET')
+        ax[i, j].grid('on')
+
+    fig.tight_layout()
+    # ------------------------------------------------------------------
+
+    return diff
+
+
+# ----------------------------------------------------------------------
+
+# ======================================================================
+# CL CALIBRATION FOR A NEW TRACER
+# ======================================================================
+
+
+def calib_tracer(outpib, outnew, xystr=None, figsize=(10, 10), fontsize=12):
+    '''calibrate CL scale for a new amyloid PET tracer,
+       given the dictionaries from processing the PiB
+       and a new tracer scans, i.e., `outpib` and `outnew`.
+
+       Other arguments:
+       xystr - dictionary of plot string XY locations.
+       figsize - tuple for controlling the figure size
+    '''
+
+    # > get PiB centiloid anchor points
+    clap = get_cl_anchors()
+
+    if xystr is None:
+        xystr = {'wc': [1.5, 1.1], 'cg': [2.0, 1.2], 'wcb': [1.5, 1.0], 'pns': [1.0, 0.75]}
+
+    # > extracting the index from names
+    pp = re.compile(r'(GE_\w*_\w*_(\d*)_NIFTI)|(Y*(\d*)\w*)')
+
+    # > calibration dictionary with SUVr and CL values
+    cal = {}
+    for rvoi in rvois:
+
+        print('------------------------------------------------------------')
+        print(rvoi_str[rvoi])
+
+        # > PiB and F18 SUVrs and PiB CL in one list/array
+        cl_suvr = []
+
+        # > initialise dictionary
+        cal[rvoi] = {'sbj': {}, 'calib': {}}
+
+        for k in outpib:
+
+            idx = pp.match(k)[2] or pp.match(k)[4]
+
+            # > get the key name in the NEW tracer dataset
+            # > (matching using the first 20 characters)
+            kf = [ky for ky in outnew if idx in ky[:20]] # idx+'_'
+            if len(kf) == 1:
+                kf = kf[0]
+            else:
+                raise ValueError('e> problem with F18 index...')
+
+            # > get the SUVr for PiB and FBB
+            suvrp = outpib[k]['suvr'][rvoi]
+            suvrf = outnew[kf]['suvr'][rvoi]
+
+            # > calculate the CL for the calibrating PiB scans
+            cl = 100 * (suvrp - clap[rvoi][0]) / (clap[rvoi][1] - clap[rvoi][0])
+
+            # > put the SUVrs for PiB and FBB into dictionary
+            cal[rvoi]['sbj'][idx] = {'cl': cl, 'suvrp': suvrp, 'suvrf': suvrf}
+
+            # > ... and a list for a future array
+            cl_suvr.append([cl, cal[rvoi]['sbj'][idx]['suvrp'], cal[rvoi]['sbj'][idx]['suvrf']])
+
+            print(f'refvoi={rvoi}, indx> {idx}, cl={cl:.3f},'
+                  f' suvr_pib={suvrp:.3f}, suvr_fbb={suvrf:.3f}')
+
+        # --------------------------------------------------------------
+        # > find the linear relationship for FBB_SUVr = m*PiB_SUVr + b
+        # > calculate the tracer `m_std` and `b_std`
+        # > (Eq. 2.2.3.1a in Klunk et al. 2015)
+
+        cl_suvr = np.array(cl_suvr)
+        m_std, b_std, r, p, stderr = linregress(cl_suvr[:, 1], cl_suvr[:, 2])
+        r2 = r**2
+
+        suvr_pib_calc = (cl_suvr[:, 2] - b_std) / m_std
+
+        cl_std_fbb = 100 * (suvr_pib_calc - clap[rvoi][0]) / (clap[rvoi][1] - clap[rvoi][0])
+
+        cal[rvoi]['calib'] = {
+            'm_std': m_std, 'b_std': b_std, 'r2': r2, 'p': p, 'stderr': stderr, 'cl_suvr': cl_suvr,
+            'suvr_pib_calc': suvr_pib_calc, 'cl_std_fbb': cl_std_fbb}
+        # --------------------------------------------------------------
+
+    # ------------------------------------------------------------------
+    # VISUALISATION
+
+    fig, ax = plt.subplots(2, 2, figsize=figsize)
+    fig2, ax2 = plt.subplots(2, 2, figsize=figsize)
+
+    for a, rvoi in enumerate(rvois):
+        i = a // 2
+        j = a % 2
+
+        ax[i, j].scatter(cal[rvoi]['calib']['cl_suvr'][:, 1], cal[rvoi]['calib']['cl_suvr'][:, 2],
+                         c='black')
+        identity_line(ax=ax[i, j], ls='--', c='b')
+        m_std = cal[rvoi]['calib']['m_std']
+        b_std = cal[rvoi]['calib']['b_std']
+        r2 = cal[rvoi]['calib']['r2']
+        ax[i, j].text(xystr[rvoi][0], xystr[rvoi][1], f'$y = {m_std:.4f}x + {b_std:.4f}$',
+                      fontsize=fontsize)
+        ax[i, j].text(xystr[rvoi][0], xystr[rvoi][1] - 0.1, f'$R^2={r2:.4f}$', fontsize=fontsize)
+        ax[i, j].set_title(rvoi_str[rvoi], fontweight='bold')
+        ax[i, j].set_xlabel(r'$^\mathrm{PiB}$SUVr$_\mathrm{IND}$', fontsize=fontsize)
+        ax[i, j].set_ylabel(r'$^\mathrm{NEW}$SUVr$_\mathrm{IND}$', fontsize=fontsize)
+        ax[i, j].grid('on')
+
+        low_x, high_x = ax[i, j].get_xlim()
+        x = np.linspace(low_x, high_x)
+        y = m_std*x + b_std
+        ax[i, j].plot(x, y, 'g')
+
+        ax2[i, j].scatter(cal[rvoi]['calib']['cl_suvr'][:, 0], cal[rvoi]['calib']['cl_std_fbb'],
+                          c='black')
+        identity_line(ax=ax2[i, j], ls='--', c='b')
+        ax2[i, j].set_xlabel('CL$^{**}$', fontsize=fontsize)
+        ax2[i, j].set_ylabel(r'$^\mathrm{NEW}$CL$_\mathrm{Std}$', fontsize=fontsize)
+        ax2[i, j].grid('on')
+        ax2[i, j].set_title(rvoi_str[rvoi])
+        # > add the same text about the equation used for converting to PiB SUVrs
+        ax2[i, j].text(
+            40, 10, f'$y_{{\\mathrm{{SUVr}}}}^{{\\mathrm{{NEW}}}}'
+            f' = {m_std:.4f}x_{{\\mathrm{{SUVr}}}}^{{\\mathrm{{PiB}}}} + {b_std:.4f}$',
+            fontsize=fontsize)
+        ax2[i, j].text(40, -5, f'$R^2={r2:.4f}$', fontsize=fontsize)
+
+    fig.tight_layout()
+    fig2.tight_layout()
+    # ------------------------------------------------------------------
+
+    return cal
+
+
+# ----------------------------------------------------------------------
+def save_cl_anchors(diff, outpath=None):
+    """save the CL anchor points for each reference VOI."""
+
+    CLA = {}
+    txt = '# centiloid anchor points for different reference VOIs\n'
+    for rvoi in rvois:
+        suvr_0 = diff['yc'][rvoi]['mean']
+        suvr_100 = diff['ad'][rvoi]['mean']
+        CLA[rvoi] = (suvr_0, suvr_100)
+        txt += f'cla_{rvoi} = ({suvr_0}, {suvr_100})\n'
+    print(txt)
+
+    if outpath is None:
+        outpath = cl_anchor_fldr
+
+    outpath = Path(outpath) / 'CL_PiB_anchors.pkl'
+    nimpa.create_dir(outpath.parent)
+    with open(outpath, 'wb') as f:
+        pickle.dump(CLA, f)
+
+    CLA['path'] = outpath
+    return CLA
+
+
+def get_cl_anchors(path=None):
+    """load the centiloid anchor points"""
+    if path is None:
+        path = cl_anchor_fldr
+    path = Path(path) / 'CL_PiB_anchors.pkl'
+    with open(path, 'rb') as f:
+        return pickle.load(f)
+
+
+# ----------------------------------------------------------------------
+
+
+# ----------------------------------------------------------------------
+def save_suvr2pib(cal, tracer, outpath=None):
+    ''' save the linear transformation parameters, `m_Std` and `b_Std`;
+        the parameters are taken from the bigger calibration dictionary `cal`
+        used here as input.
+
+        The compact dictionary consists the two parameters for all reference regions
+        and for the specific F-18 tracer being used.
+
+        `tracer` can be:
+            - 'fbb' for [18F]florbetaben
+            - 'flute' for [18F]flutemetamol
+            - 'fbp' for [18F]florbetapir
+    '''
+
+    if tracer not in ('fbp', 'fbb', 'flute'):
+        raise ValueError('tracer is unrecognised or not given!')
+
+    # conversion dictionary
+    CNV = {}
+    txt = '# Tracer SUVr conversion to PiB SUVr\n'
+    for rvoi in rvois:
+        mstd = cal[rvoi]['calib']['m_std']
+        bstd = cal[rvoi]['calib']['b_std']
+        CNV[rvoi] = {'m_std': mstd, 'b_std': bstd}
+        txt += f'{rvoi}: m_std, b_std = ({mstd}, {bstd})\n'
+    print(txt)
+
+    if outpath is None:
+        outpath = cl_anchor_fldr
+    outpath = Path(outpath)
+
+    nimpa.create_dir(outpath)
+    pth = outpath / f'suvr_{tracer}_to_suvr_pib__transform.pkl'
+    with open(pth, 'wb') as f:
+        pickle.dump(CNV, f)
+
+    CNV['path'] = pth
+    return CNV
+
+
+def get_suvr2pib(tracer, path=None):
+    """
+    load the linear transformation parameters from a tracer SUVr
+    to PiB SUVr; it's used for converting to CL scale any F-18
+    tracers.
+    """
+    if tracer not in ['fbp', 'fbb', 'flute']:
+        raise ValueError(f'tracer ({tracer}) unrecognised')
+
+    if path is None:
+        path = cl_anchor_fldr
+    pth = Path(path) / f'suvr_{tracer}_to_suvr_pib__transform.pkl'
+
+    with open(pth, 'rb') as f:
+        return pickle.load(f)
diff --git a/amypet/backend_centiloid.py b/amypet/backend_centiloid.py
new file mode 100644
index 00000000..dc853949
--- /dev/null
+++ b/amypet/backend_centiloid.py
@@ -0,0 +1,461 @@
+"""Centiloid pipeline
+
+Usage:
+  centiloid [options] <fpets> <fmris>
+
+Arguments:
+  <fpets>  : PET NIfTI directory [default: DirChooser]
+  <fmris>  : MRI NIfTI directory [default: DirChooser]
+
+Options:
+  --outpath FILE  : Output directory
+  --visual  : whether to plot
+"""
+
+__author__ = ("Pawel J Markiewicz", "Casper da Costa-Luis")
+__copyright__ = "Copyright 2022"
+
+import logging
+import os
+import pickle
+from pathlib import Path
+from typing import Optional
+
+import matplotlib.pyplot as plt
+import numpy as np
+import spm12
+from miutil.fdio import hasext, nsort
+from niftypet import nimpa
+
+from .aux import cl_anchor_fldr, cl_masks_fldr
+
+log = logging.getLogger(__name__)
+
+
+# ----------------------------------------------------------------------
+def load_masks(mskpath, voxsz: int = 2):
+    ''' Load the Centiloid PET masks for calculating
+        the SUVr to then convert it to Centiloid.
+
+        Return the paths to the masks and the masks themselves
+        Options:
+        - voxsz: voxel size used in the normalised PET/MR images
+                 and the atlas mask.  Default is 2 mm.
+    '''
+
+    voxsz = int(voxsz)
+
+    if voxsz not in [1, 2]:
+        raise ValueError('Incorrect voxel size - only 1 and 2 are accepted.')
+
+    log.info('loading CL masks...')
+    fmasks = {
+        'cg': mskpath / f'voi_CerebGry_{voxsz}mm.nii', 'wc': mskpath / f'voi_WhlCbl_{voxsz}mm.nii',
+        'wcb': mskpath / f'voi_WhlCblBrnStm_{voxsz}mm.nii',
+        'pns': mskpath / f'voi_Pons_{voxsz}mm.nii', 'ctx': mskpath / f'voi_ctx_{voxsz}mm.nii'}
+    masks = {fmsk: nimpa.getnii(fmasks[fmsk]) for fmsk in fmasks}
+
+    return fmasks, masks
+
+
+def run(fpets, fmris, tracer='pib', flip_pet=None, bias_corr=True, voxsz: int = 2, outpath=None,
+        visual=False, climage=True, used_saved=False, cl_anchor_path: Optional[Path] = None):
+    """
+    Process centiloid (CL) using input file lists for PET and MRI
+    images, `fpets` and `fmris` (must be in NIfTI format).
+    Args:
+      outpath: path to the output folder
+      bias_corr: bias filed correction for the MR image (True/False)
+      tracer: specifies what tracer is being used and so the right
+              transformation is used; by default it is PiB.  Currently
+              [18F]flutemetamol, 'flute', and [18F]florbetaben, 'fbb'
+              are supported.
+              IMPORTANT: when calibrating a new tracer, ensure that
+              `tracer`='new'.
+      voxsz: voxel size for SPM normalisation writing function
+             (output MR and PET images will have this voxel size).
+      flip_pet: a list of flips (3D tuples) which flip any dimension
+               of the 3D PET image (z,y,x); the list has to have the
+               same length as the lists of `fpets` and `fmris`
+      used_saved: if True, looks for already saved normalised PET
+                images and loads them to avoid processing time.
+      visual: SPM-based progress visualisations of image registration or
+              or image normalisation.
+      climage: outputs the CL-converted PET image in the MNI space
+      cl_anchor_path: The path where optional CL anchor dictionary is
+                saved.
+    """
+
+    # supported F-18 tracers
+    f18_tracers = ['fbp', 'fbb', 'flute']
+
+    spm_path = Path(spm12.utils.spm_dir())
+
+    out = {}                                           # output dictionary
+    tmpl_avg = spm_path / 'canonical' / 'avg152T1.nii' # template path
+
+    if isinstance(fpets, (str, Path)) and isinstance(fmris, (str, Path)):
+        # when single PET and MR files are given
+        if os.path.isfile(fpets) and os.path.isfile(fmris):
+            pet_mr_list = [[fpets], [fmris]]
+        # when folder paths are given for PET and MRI files
+        # the content is sorted for both PET and MRI and then matched
+        # according to the order
+        elif os.path.isdir(fpets) and os.path.isdir(fmris):
+            fp = nsort(
+                os.path.join(fpets, f) for f in os.listdir(fpets) if hasext(f, ('nii', 'nii.gz')))
+            fm = nsort(
+                os.path.join(fmris, f) for f in os.listdir(fmris) if hasext(f, ('nii', 'nii.gz')))
+            pet_mr_list = [fp, fm]
+        else:
+            raise ValueError('unrecognised or unmatched input for PET and MRI image data')
+    elif isinstance(fpets, list) and isinstance(fmris, list):
+        if len(fpets) != len(fmris):
+            raise ValueError('the number of PET and MRI files have to match!')
+        # check if all files exists
+        if not all(os.path.isfile(f) and hasext(f, ('nii', 'nii.gz')) for f in fpets) or not all(
+                os.path.isfile(f) and hasext(f, ('nii', 'nii.gz')) for f in fmris):
+            raise ValueError('the number of paired files do not match')
+        pet_mr_list = [fpets, fmris]
+    else:
+        raise ValueError('unrecognised input image data')
+
+    # -------------------------------------------------------------
+    if flip_pet is not None:
+        if isinstance(flip_pet, tuple) and len(pet_mr_list[0]) == 1:
+            flips = [flip_pet]
+        elif isinstance(flip_pet, list) and len(flip_pet) == len(pet_mr_list[0]):
+            flips = flip_pet
+        else:
+            log.warning('the flip definition is not compatible with the list of PET images')
+    else:
+        flips = [None] * len(pet_mr_list[0])
+    # -------------------------------------------------------------
+
+    # -------------------------------------------------------------
+    # > get the CL masks
+    fmasks, masks = load_masks(cl_masks_fldr, voxsz=voxsz)
+    # -------------------------------------------------------------
+
+    log.info('iterate through all the input data...')
+
+    fi = -1
+    for fpet, fmri in zip(*pet_mr_list):
+
+        fi += 1
+        # > make the files Path objects
+        fpet, fmri = Path(fpet), Path(fmri)
+        opth = Path(fpet.parent if outpath is None else outpath)
+
+        # > create output dictionary and folder specific to the PET
+        # > file which will be CL quantified; name of the folder
+        # > and dictionary output is based on PET name
+        onm = fpet.name.rsplit('.nii', 1)[0]
+        spth = opth / onm
+        out[onm] = {'opth': spth, 'fpet': fpet, 'fmri': fmri}
+
+        # >output path for centre of mass alignment and registration
+        opthc = spth / 'centre-of-mass'
+        opthr = spth / 'registration'
+        opthn = spth / 'normalisation'
+        optho = spth / 'normalised'
+        opths = spth / 'suvr'
+        opthi = spth / 'cl-image'
+
+        # > find if the normalised PET is already there
+        if optho.is_dir():
+            fnpets = [f for f in optho.iterdir() if fpet.name.split('.nii')[0] in f.name]
+        else:
+            fnpets = []
+
+        if used_saved and len(fnpets) == 1:
+            log.info(f'subject {onm}: loading already normalised PET image...')
+
+        else:
+
+            # run bias field correction unless cancelled
+            if bias_corr:
+                log.info(f'subject {onm}: running MR bias field correction')
+                out[onm]['n4'] = nimpa.bias_field_correction(fmri, executable='sitk', outpath=spth)
+                fmri = out[onm]['n4']['fim']
+
+            log.info(f'subject {onm}: centre of mass correction')
+            # > check if flipping the PET is requested
+            if flips[fi] is not None and any(flips[fi]):
+                flip = flips[fi]
+            else:
+                flip = None
+
+            # > modify for the centre of mass being at O(0,0,0)
+            out[onm]['petc'] = petc = nimpa.centre_mass_corr(fpet, flip=flip, outpath=opthc)
+            out[onm]['mric'] = mric = nimpa.centre_mass_corr(fmri, outpath=opthc)
+
+            log.info(f'subject {onm}: MR registration to MNI space')
+            out[onm]['reg1'] = reg1 = spm12.coreg_spm(
+                tmpl_avg,
+                mric['fim'],
+                fwhm_ref=0,
+                fwhm_flo=3,
+                outpath=opthr,
+                fname_aff="",
+                fcomment="",
+                pickname="ref",
+                costfun="nmi",
+                graphics=1,
+                visual=int(visual),
+                del_uncmpr=True,
+                save_arr=True,
+                save_txt=True,
+                modify_nii=True,
+            )
+
+            log.info(f'subject {onm}: PET -> MR registration')
+            out[onm]['reg2'] = reg2 = spm12.coreg_spm(
+                reg1['freg'],
+                petc['fim'],
+                fwhm_ref=3,
+                fwhm_flo=6,
+                outpath=opthr,
+                fname_aff="",
+                fcomment='_mr-reg',
+                pickname="ref",
+                costfun="nmi",
+                graphics=1,
+                visual=int(visual),
+                del_uncmpr=True,
+                save_arr=True,
+                save_txt=True,
+                modify_nii=True,
+            )
+
+            log.info(f'subject {onm}: MR normalisation/segmentation...')
+            out[onm]['norm'] = norm = spm12.seg_spm(reg1['freg'], spm_path, outpath=opthn,
+                                                    store_nat_gm=True, store_nat_wm=False,
+                                                    store_nat_csf=True, store_fwd=True,
+                                                    store_inv=True, visual=visual)
+            # > normalise
+            list4norm = [reg1['freg'] + ',1', reg2['freg'] + ',1']
+            out[onm]['fnorm'] = spm12.normw_spm(norm['fordef'], list4norm, voxsz=float(voxsz),
+                                                outpath=optho)
+
+            log.info(f'subject {onm}: load normalised PET image...')
+            fnpets = [
+                f for f in optho.iterdir()
+                if fpet.name.split('.nii')[0] in f.name and 'n4bias' not in f.name.lower()]
+            # and 'mr' not in f.name.lower()
+
+            if len(fnpets) == 0:
+                raise ValueError('could not find normalised PET image files')
+            elif len(fnpets) > 1:
+                raise ValueError('too many potential normalised PET image files found')
+
+        npet_dct = nimpa.getnii(fnpets[0], output='all')
+        npet = npet_dct['im']
+        npet[np.isnan(npet)] = 0 # get rid of NaNs
+
+        # npet[npet<0] = 0
+
+        # > extract mean values and SUVr
+        out[onm]['avgvoi'] = avgvoi = {fmsk: np.mean(npet[masks[fmsk] > 0]) for fmsk in fmasks}
+        out[onm]['suvr'] = suvr = {
+            fmsk: avgvoi['ctx'] / avgvoi[fmsk]
+            for fmsk in fmasks if fmsk != 'ctx'}
+
+        # **************************************************************
+        # C E N T I L O I D   S C A L I N G
+        # **************************************************************
+        # ---------------------------------
+        # > path to anchor point dictionary
+        if cl_anchor_path is None:
+            cl_fldr = cl_anchor_fldr
+        else:
+            cl_fldr = Path(cl_anchor_path)
+        assert cl_fldr.is_dir()
+        # ---------------------------------
+
+        # ---------------------------------
+        # > centiloid transformation for PiB
+        pth = cl_fldr / 'CL_PiB_anchors.pkl'
+
+        if not os.path.isfile(pth):
+            tracer = 'new'
+            log.warning('Could not find the PiB CL anchor point definitions/tables')
+        else:
+            log.info(f'using the following CL anchor table:\n   {pth}')
+
+        with open(pth, 'rb') as f:
+            CLA = pickle.load(f)
+        # ---------------------------------
+
+        # ---------------------------------
+        # > centiloid transformation for PiB
+        # > check if SUVr transformation is needed for F-18 tracers
+        if tracer in f18_tracers:
+            pth = cl_fldr / f'suvr_{tracer}_to_suvr_pib__transform.pkl'
+
+            if not os.path.isfile(pth):
+                log.warning(
+                    'The conversion dictionary/table for the specified tracer is not found')
+            else:
+                with open(pth, 'rb') as f:
+                    CNV = pickle.load(f)
+                print('loaded SUVr transformations for tracer:', tracer)
+
+                # > save the PiB converted SUVrs
+                out[onm]['suvr_pib_calc'] = suvr_pib_calc = {
+                    fmsk: (suvr[fmsk] - CNV[fmsk]['b_std']) / CNV[fmsk]['m_std']
+                    for fmsk in fmasks if fmsk != 'ctx'}
+
+                # > save the linear transformation parameters
+                out[onm]['suvr_pib_calc_transf'] = {
+                    fmsk: (CNV[fmsk]['m_std'], CNV[fmsk]['b_std'])
+                    for fmsk in fmasks if fmsk != 'ctx'}
+
+            # > used now the new PiB converted SUVrs
+            suvr = suvr_pib_calc
+        # ---------------------------------
+
+        if tracer != 'new':
+            out[onm]['cl'] = cl = {
+                fmsk: 100 * (suvr[fmsk] - CLA[fmsk][0]) / (CLA[fmsk][1] - CLA[fmsk][0])
+                for fmsk in fmasks if fmsk != 'ctx'}
+
+        # **************************************************************
+
+        # -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
+        # > output the CL-converted PET image in the MNI space
+        if climage and tracer != 'new':
+
+            for refvoi in fmasks:
+                if refvoi == 'ctx':
+                    continue
+
+                refavg = out[onm]['avgvoi'][refvoi]
+
+                # > obtain the SUVr image for the given reference VOI
+                npet_suvr = npet / refavg
+
+                # > convert to PiB scale if it is an F18 tracer
+                if tracer in f18_tracers:
+                    npet_suvr = (npet_suvr - CNV[refvoi]['b_std']) / CNV[refvoi]['m_std']
+
+                # > convert the (PiB) SUVr image to CL scale
+                npet_cl = 100 * (npet_suvr - CLA[refvoi][0]) / (CLA[refvoi][1] - CLA[refvoi][0])
+
+                # > get the CL global value by applying the CTX mask
+                cl_ = np.mean(npet_cl[masks['ctx'].astype(bool)])
+
+                cl_refvoi = cl[refvoi]
+                if tracer != 'new' and abs(cl_ - cl_refvoi) < 0.25:
+                    # > save the CL-converted file
+                    nimpa.create_dir(opthi)
+                    fout = opthi / ('CL_image_ref-' + refvoi + fnpets[0].name.split('.nii')[0] +
+                                    '.nii.gz')
+                    nimpa.array2nii(
+                        npet_cl, npet_dct['affine'], fout,
+                        trnsp=(npet_dct['transpose'].index(0), npet_dct['transpose'].index(1),
+                               npet_dct['transpose'].index(2)), flip=npet_dct['flip'])
+
+                elif tracer != 'new' and abs(cl_ - cl_refvoi) > 0.25:
+                    log.warning('The CL of CL-converted image is different to the calculated CL'
+                                f' (CL_img={cl_:.4f} vs CL={cl_refvoi:.4f}).')
+                    log.warning('The CL image has not been generated!')
+                else:
+                    log.warning('The CL image has not been generated due to new tracer being used')
+        # -*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
+
+        # -------------------------------------------------------------
+        # VISUALISATION
+        # pick masks for visualisation
+        msk = 'ctx'
+        mskr = 'wc' # 'pons'#
+
+        showpet = nimpa.imsmooth(npet.astype(np.float32), voxsize=npet_dct['voxsize'], fwhm=3.)
+
+        nimpa.create_dir(opths)
+
+        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+        # transaxial tiling for QC
+        # choose the shape of the mosaic/tiled image
+        izs = np.array([[30, 40, 50], [60, 70, 80]])
+
+        # get the shape of the mosaic
+        shp = izs.shape
+
+        # initialise mosaic for PET and masks
+        mscp_t = np.zeros((shp[0], shp[1]) + showpet[0, ...].shape, dtype=np.float32)
+        mscm_t = mscp_t.copy()
+
+        # fill in the images
+        for i in range(shp[0]):
+            for j in range(shp[1]):
+                mscp_t[i, j, ...] = showpet[izs[i, j], ...]
+                mscm_t[i, j, ...] = masks[msk][izs[i, j], ...] + masks[mskr][izs[i, j], ...]
+
+        # reshape for the final touch
+        mscp_t = mscp_t.swapaxes(1, 2)
+        mscm_t = mscm_t.swapaxes(1, 2)
+        mscp_t = mscp_t.reshape(shp[0] * showpet.shape[1], shp[1] * showpet.shape[2])
+        mscm_t = mscm_t.reshape(shp[0] * showpet.shape[1], shp[1] * showpet.shape[2])
+
+        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+        # Sagittal tiling for QC
+
+        # choose the shape of the mosaic/tiled image
+        ixs = np.array([[20, 30, 40], [50, 60, 70]])
+
+        # get the shape of the mosaic
+        shp = ixs.shape
+
+        # initialise mosaic for PET and masks
+        mscp_s = np.zeros((shp[0], shp[1]) + showpet[..., 0].shape, dtype=np.float32)
+        mscm_s = mscp_s.copy()
+
+        # fill in the images
+        for i in range(shp[0]):
+            for j in range(shp[1]):
+                mscp_s[i, j, ...] = showpet[..., ixs[i, j]]
+                mscm_s[i, j, ...] = masks[msk][..., ixs[i, j]] + masks[mskr][..., ixs[i, j]]
+
+        # reshape for the final touch
+        mscp_s = mscp_s.swapaxes(1, 2)
+        mscm_s = mscm_s.swapaxes(1, 2)
+        mscp_s = mscp_s.reshape(shp[0] * showpet.shape[0], shp[1] * showpet.shape[1])
+        mscm_s = mscm_s.reshape(shp[0] * showpet.shape[0], shp[1] * showpet.shape[1])
+
+        # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+        fig, ax = plt.subplots(2, 1, figsize=(9, 12))
+
+        thrsh = 0.9 * showpet.max()
+
+        ax[0].imshow(mscp_t, cmap='magma', vmax=thrsh)
+        ax[0].imshow(mscm_t, cmap='gray_r', alpha=0.25)
+        ax[0].set_axis_off()
+        ax[0].set_title(f'{onm}: transaxial centiloid sampling')
+
+        ax[1].imshow(mscp_s, cmap='magma', vmax=thrsh)
+        ax[1].imshow(mscm_s, cmap='gray_r', alpha=0.25)
+        ax[1].set_axis_off()
+        ax[1].set_title(f'{onm} sagittal centiloid sampling')
+
+        suvrstr = ",   ".join([
+            f"PiB transformed: $SUVR_{{WC}}=${suvr['wc']:.3f}", f"$SUVR_{{GMC}}=${suvr['cg']:.3f}",
+            f"$SUVR_{{CBS}}=${suvr['wcb']:.3f}", f"$SUVR_{{PNS}}=${suvr['pns']:.3f}"])
+        ax[1].text(0, 190, suvrstr, fontsize=12)
+
+        if tracer != 'new':
+            clstr = ",   ".join([
+                f"$CL_{{WC}}=${cl['wc']:.1f}", f"$CL_{{GMC}}=${cl['cg']:.1f}",
+                f"$CL_{{CBS}}=${cl['wcb']:.1f}", f"$CL_{{PNS}}=${cl['pns']:.1f}"])
+            ax[1].text(0, 205, clstr, fontsize=12)
+
+        plt.tight_layout()
+
+        fqcpng = opths / (onm+'_CL-SUVr_mask_PET_sampling.png')
+        plt.savefig(fqcpng, dpi=150, facecolor='auto', edgecolor='auto')
+
+        plt.close('all')
+
+        out[onm]['fqc'] = fqcpng
+
+    return out
diff --git a/amypet/centiloid.py b/amypet/centiloid.py
new file mode 100644
index 00000000..d5c6a190
--- /dev/null
+++ b/amypet/centiloid.py
@@ -0,0 +1,31 @@
+"""Centiloid pipeline
+
+Usage:
+  centiloid [options] <inpath>
+
+Arguments:
+  <inpath>  : PET & MRI NIfTI directory [default: DirChooser]
+
+Options:
+  --outpath DIR  : Output directory
+  --tracer TYPE  : PET Tracer; choices: {[default: pib],flt,fbp,fbb}
+  --start TIME  : SUVr start time [default: None:float]
+  --end TIME  : SUVr end time [default: None:float]
+  --dynamic-analysis  : whether to do dynamic analysis
+"""
+__author__ = "Casper da Costa-Luis"
+__copyright__ = "Copyright 2022"
+import logging
+
+from .backend_centiloid import run as centiloid_run
+
+log = logging.getLogger(__name__)
+TRACERS = ('pib', 'flt', 'fbp', 'fbb')
+
+
+def run(inpath, tracer='pib', start=None, end=None, dynamic_analysis=False, outpath=None):
+    """Just a stub"""
+    # probably want to use `centiloid_run()`?
+    return {
+        'inpath': inpath, 'tracer': tracer, 'start': start, 'end': end,
+        'dynamic': dynamic_analysis, 'outpath': outpath}
diff --git a/amypet/centiloid/__init__.py b/amypet/centiloid_old/__init__.py
similarity index 100%
rename from amypet/centiloid/__init__.py
rename to amypet/centiloid_old/__init__.py
diff --git a/amypet/centiloid/f_1CorregisterEstimate.m b/amypet/centiloid_old/f_1CorregisterEstimate.m
similarity index 100%
rename from amypet/centiloid/f_1CorregisterEstimate.m
rename to amypet/centiloid_old/f_1CorregisterEstimate.m
diff --git a/amypet/centiloid/f_2CorregisterEstimate.m b/amypet/centiloid_old/f_2CorregisterEstimate.m
similarity index 100%
rename from amypet/centiloid/f_2CorregisterEstimate.m
rename to amypet/centiloid_old/f_2CorregisterEstimate.m
diff --git a/amypet/centiloid/f_3Segment.m b/amypet/centiloid_old/f_3Segment.m
similarity index 100%
rename from amypet/centiloid/f_3Segment.m
rename to amypet/centiloid_old/f_3Segment.m
diff --git a/amypet/centiloid/f_4Normalise.m b/amypet/centiloid_old/f_4Normalise.m
similarity index 100%
rename from amypet/centiloid/f_4Normalise.m
rename to amypet/centiloid_old/f_4Normalise.m
diff --git a/amypet/centiloid/f_Quant_centiloid.m b/amypet/centiloid_old/f_Quant_centiloid.m
similarity index 100%
rename from amypet/centiloid/f_Quant_centiloid.m
rename to amypet/centiloid_old/f_Quant_centiloid.m
diff --git a/amypet/centiloid/f_acpcReorientation.m b/amypet/centiloid_old/f_acpcReorientation.m
similarity index 100%
rename from amypet/centiloid/f_acpcReorientation.m
rename to amypet/centiloid_old/f_acpcReorientation.m
diff --git a/amypet/centiloid/f_noNaN.m b/amypet/centiloid_old/f_noNaN.m
similarity index 100%
rename from amypet/centiloid/f_noNaN.m
rename to amypet/centiloid_old/f_noNaN.m
diff --git a/amypet/data/CL_masks/voi_CerebGry_1mm.nii b/amypet/data/CL_masks/voi_CerebGry_1mm.nii
new file mode 100644
index 00000000..07b53a32
Binary files /dev/null and b/amypet/data/CL_masks/voi_CerebGry_1mm.nii differ
diff --git a/amypet/data/CL_masks/voi_CerebGry_2mm.nii b/amypet/data/CL_masks/voi_CerebGry_2mm.nii
new file mode 100644
index 00000000..07706e52
Binary files /dev/null and b/amypet/data/CL_masks/voi_CerebGry_2mm.nii differ
diff --git a/amypet/data/CL_masks/voi_Pons_1mm.nii b/amypet/data/CL_masks/voi_Pons_1mm.nii
new file mode 100644
index 00000000..5a17bbb6
Binary files /dev/null and b/amypet/data/CL_masks/voi_Pons_1mm.nii differ
diff --git a/amypet/data/CL_masks/voi_Pons_2mm.nii b/amypet/data/CL_masks/voi_Pons_2mm.nii
new file mode 100644
index 00000000..af039692
Binary files /dev/null and b/amypet/data/CL_masks/voi_Pons_2mm.nii differ
diff --git a/amypet/data/CL_masks/voi_WhlCblBrnStm_1mm.nii b/amypet/data/CL_masks/voi_WhlCblBrnStm_1mm.nii
new file mode 100644
index 00000000..1bc8cabe
Binary files /dev/null and b/amypet/data/CL_masks/voi_WhlCblBrnStm_1mm.nii differ
diff --git a/amypet/data/CL_masks/voi_WhlCblBrnStm_2mm.nii b/amypet/data/CL_masks/voi_WhlCblBrnStm_2mm.nii
new file mode 100644
index 00000000..50a45b02
Binary files /dev/null and b/amypet/data/CL_masks/voi_WhlCblBrnStm_2mm.nii differ
diff --git a/amypet/data/CL_masks/voi_WhlCbl_1mm.nii b/amypet/data/CL_masks/voi_WhlCbl_1mm.nii
new file mode 100644
index 00000000..b5625e1c
Binary files /dev/null and b/amypet/data/CL_masks/voi_WhlCbl_1mm.nii differ
diff --git a/amypet/data/CL_masks/voi_WhlCbl_2mm.nii b/amypet/data/CL_masks/voi_WhlCbl_2mm.nii
new file mode 100644
index 00000000..cb44aa44
Binary files /dev/null and b/amypet/data/CL_masks/voi_WhlCbl_2mm.nii differ
diff --git a/amypet/data/CL_masks/voi_ctx_1mm.nii b/amypet/data/CL_masks/voi_ctx_1mm.nii
new file mode 100644
index 00000000..41f3ef3f
Binary files /dev/null and b/amypet/data/CL_masks/voi_ctx_1mm.nii differ
diff --git a/amypet/data/CL_masks/voi_ctx_2mm.nii b/amypet/data/CL_masks/voi_ctx_2mm.nii
new file mode 100644
index 00000000..7dea9fda
Binary files /dev/null and b/amypet/data/CL_masks/voi_ctx_2mm.nii differ
diff --git a/amypet/data/cl_anchor_tables/CL_PiB_anchors.pkl b/amypet/data/cl_anchor_tables/CL_PiB_anchors.pkl
new file mode 100644
index 00000000..5f969d2e
Binary files /dev/null and b/amypet/data/cl_anchor_tables/CL_PiB_anchors.pkl differ
diff --git a/amypet/data/cl_anchor_tables/suvr_fbb_to_suvr_pib__transform.pkl b/amypet/data/cl_anchor_tables/suvr_fbb_to_suvr_pib__transform.pkl
new file mode 100644
index 00000000..c65d0596
Binary files /dev/null and b/amypet/data/cl_anchor_tables/suvr_fbb_to_suvr_pib__transform.pkl differ
diff --git a/amypet/data/cl_anchor_tables/suvr_fbp_to_suvr_pib__transform.pkl b/amypet/data/cl_anchor_tables/suvr_fbp_to_suvr_pib__transform.pkl
new file mode 100644
index 00000000..c4b9b84c
Binary files /dev/null and b/amypet/data/cl_anchor_tables/suvr_fbp_to_suvr_pib__transform.pkl differ
diff --git a/amypet/data/cl_anchor_tables/suvr_flute_to_suvr_pib__transform.pkl b/amypet/data/cl_anchor_tables/suvr_flute_to_suvr_pib__transform.pkl
new file mode 100644
index 00000000..e610edd6
Binary files /dev/null and b/amypet/data/cl_anchor_tables/suvr_flute_to_suvr_pib__transform.pkl differ
diff --git a/amypet/gui.py b/amypet/gui.py
index 7a593311..9050cea9 100755
--- a/amypet/gui.py
+++ b/amypet/gui.py
@@ -257,8 +257,12 @@ def argparser(prog, description=None, epilog=None, formatter_class=None):
     Func(imtrimup.run, imtrimup.__doc__, version=niftypad.__version__, python_deps=["nimpa"],
          argparser=argparser)
 
-    Func(centiloid.run, centiloid.__doc__, version=niftypad.__version__,
-         python_deps=["miutil[nii]", "setuptools", "spm12", "tqdm"], argparser=argparser)
+    centiloid_func = Func(centiloid.run, centiloid.__doc__, version=niftypad.__version__,
+                          python_deps=["miutil[nii]", "nimpa", "setuptools", "spm12",
+                                       "tqdm"], argparser=argparser)
+    opts = centiloid_func.parser._get_optional_actions()
+    tracer = next(i for i in opts if i.dest == 'tracer')
+    tracer.choices = centiloid.TRACERS
 
     kinetic_model = Func(
         niftypad.api.kinetic_model, """\
diff --git a/amypet/preproc.py b/amypet/preproc.py
new file mode 100644
index 00000000..abf361c9
--- /dev/null
+++ b/amypet/preproc.py
@@ -0,0 +1,494 @@
+'''
+Preprocessing tools for AmyPET core processes
+'''
+
+__author__ = "Pawel Markiewicz"
+__copyright__ = "Copyright 2022"
+
+import logging as log
+import os
+import shutil
+from datetime import datetime
+from itertools import combinations
+from pathlib import Path
+from subprocess import run
+
+import dcm2niix
+import numpy as np
+import spm12
+from niftypet import nimpa
+
+from .aux import get_atlas
+from .suvr_tools import r_trimup
+
+log.basicConfig(level=log.WARNING, format=nimpa.LOG_FORMAT)
+
+# ------------------------------------------------
+# DEFINITIONS:
+# TODO: move these to a separate file, e.g., `defs.py`
+
+# > SUVr time window post injection and duration
+suvr_twindow = {
+    'flute': [90 * 60, 110 * 60, 1200], 'fbb': [90 * 60, 110 * 60, 1200],
+    'fbp': [50 * 60, 60 * 60, 600]}
+margin = 0.1
+
+# tracer names
+tracer_names = {
+    'flute': ['flt', 'flut', 'flute', 'flutemetamol'], 'fbb': ['fbb', 'florbetaben'],
+    'fbp': ['fbp', 'florbetapir']}
+
+# > break time for coffee break protocol (target)
+break_time = 1800
+
+# > time margin for the 1st coffee break acquisition
+breakdyn_t = (1200, 2400)
+
+# > minimum time for the full dynamic acquisition
+fulldyn_time = 3600
+
+# ------------------------------------------------
+
+
+# =====================================================================
+def explore_input(
+    input_fldr,
+    tracer=None,
+    suvr_win_def=None,
+    outpath=None,
+):
+    '''
+    Process the input folder of amyloid PET DICOM data.
+    The folder can contain two subfolders for a coffee break protocol including
+    early dynamic followed by a static scan.
+    The folder can also contain static or dynamic DICOM files.
+    Those files can also be within a subfolder.
+
+    Return the dictionary of (1) the list of dictionaries for each DICOM folder
+    (2) list of descriptions for each DICOM folder for classification of input
+
+    Arguments:
+    - tracer:   The name of one of the three tracers: 'flute', 'fbb', 'fbp'
+    - suvr_win_def: The definition of SUVr time frame (SUVr/CL is always calculated)
+                as a two-element list [t_start, t_stop] in seconds.  If the
+                window is not defined the function will attempt to get the
+                information from the tracer info and use the default (as
+                defined in`defs.py`)
+    - outpath:  output path where all the intermediate and final results are
+                stored.
+
+    '''
+
+    # > make the input a Path object
+    input_fldr = Path(input_fldr)
+
+    if not input_fldr.is_dir():
+        raise ValueError('Incorrect input - not a folder!')
+
+    if outpath is None:
+        amyout = input_fldr.parent / 'amypet_output'
+    else:
+        amyout = Path(outpath)
+    nimpa.create_dir(amyout)
+
+    # ================================================
+    # > first check if the folder has DICOM series
+
+    # > multiple series in folders (if any)
+    msrs = []
+    for itm in input_fldr.iterdir():
+        if itm.is_dir():
+            srs = nimpa.dcmsort(itm, grouping='a+t+d', copy_series=True, outpath=amyout)
+            if srs:
+                msrs.append(srs)
+
+    # > check files in the input folder
+    srs = nimpa.dcmsort(input_fldr, grouping='a+t+d', copy_series=True, outpath=amyout)
+    if srs:
+        msrs.append(srs)
+    # ================================================
+
+    # > initialise the list of acquisition classification
+    msrs_class = []
+
+    # > time-sorted series
+    msrs_t = []
+
+    for m in msrs:
+
+        # > for each folder do the following:
+
+        # > time sorted series according to acquisition time
+        srs_t = {k: v for k, v in sorted(m.items(), key=lambda item: item[1]['tacq'])}
+
+        msrs_t.append(srs_t)
+
+        # -----------------------------------------------
+        # > frame timings relative to the injection time -
+        #   radiopharmaceutical administration start time
+        t_frms = []
+        for k in srs_t:
+            t0 = datetime.strptime(srs_t[k]['dstudy'] + srs_t[k]['tacq'],
+                                   '%Y%m%d%H%M%S') - srs_t[k]['radio_start_time']
+            t1 = datetime.strptime(
+                srs_t[k]['dstudy'] + srs_t[k]['tacq'],
+                '%Y%m%d%H%M%S') + srs_t[k]['frm_dur'] - srs_t[k]['radio_start_time']
+            t_frms.append((t0.seconds, t1.seconds))
+
+        t_starts = [t[0] for t in t_frms]
+        t_stops = [t[1] for t in t_frms]
+
+        # > overall acquisition duration
+        acq_dur = t_frms[-1][-1] - t_frms[0][0]
+        # -----------------------------------------------
+
+        # > check if the frames qualify for static, fully dynamic or
+        # > coffee-break dynamic acquisition
+        acq_type = None
+        if t_frms[0][0] < 1:
+            if t_frms[-1][-1] > breakdyn_t[0] and t_frms[-1][-1] <= breakdyn_t[1]:
+                acq_type = 'breakdyn'
+            elif t_frms[-1][-1] >= fulldyn_time:
+                acq_type = 'fulldyn'
+        elif t_frms[0][0] > 1:
+            acq_type = 'static'
+        # -----------------------------------------------
+
+        # > classify tracer if possible and if not given
+        if tracer is None:
+            if 'tracer' in srs_t[next(iter(srs_t))]:
+                tracer_dcm = srs_t[next(iter(srs_t))]['tracer'].lower()
+                for t in tracer_names:
+                    for n in tracer_names[t]:
+                        if n in tracer_dcm:
+                            tracer = t
+
+            # > when tracer info not provided and not in DICOMs
+            if acq_type == 'static' and not tracer:
+
+                # > assuming the first of the following tracers then
+                for t in suvr_twindow:
+                    dur = suvr_twindow[t][2]
+                    if (acq_dur > dur * (1-margin)) and (t_frms[0][0] < suvr_twindow[t][0] *
+                                                         (1+margin)):
+                        tracer = t
+                        break
+        else:
+            tracer_grp = [tracer in tracer_names[t] for t in tracer_names]
+            if any(tracer_grp):
+                tracer = np.array(list(tracer_names.keys()))[tracer_grp][0]
+        # -----------------------------------------------
+
+        # > is the static acquisition covering the provided SUVr frame definition?
+        if acq_type == 'static':
+            # -----------------------------------------------
+            # > try to establish the SUVr window even if not provided
+            if suvr_win_def is None and not tracer:
+                raise ValueError(
+                    'Impossible to figure out tracer and SUVr time window - please specify them!')
+            elif suvr_win_def is None and tracer:
+                suvr_win = suvr_twindow[tracer][:2]
+            else:
+                suvr_win = suvr_win_def
+            # -----------------------------------------------
+
+            # > window margin
+            mrgn_suvr = margin * suvr_twindow[tracer][2]
+            if t_frms[0][
+                    0] < suvr_win[0] + mrgn_suvr and acq_dur > suvr_twindow[tracer][2] - mrgn_suvr:
+
+                t0_suvr = min(t_starts, key=lambda x: abs(x - suvr_win[0]))
+                t1_suvr = min(t_stops, key=lambda x: abs(x - suvr_win[1]))
+
+                frm_0 = t_starts.index(t0_suvr)
+                frm_1 = t_stops.index(t1_suvr)
+
+                msrs_class.append({
+                    'acq': [acq_type, 'suvr'], 'time': (t0_suvr, t1_suvr), 'timings': t_frms,
+                    'idxs': (frm_0, frm_1),
+                    'frms': [s for i, s in enumerate(srs_t) if i in range(frm_0, frm_1 + 1)]})
+            else:
+                log.warning('The acquisition does not cover the requested time frame!')
+
+                msrs_class.append({
+                    'acq': [acq_type], 'time': (t_starts[0], t_stops[-1]),
+                    'idxs': (0, len(t_frms) - 1), 'frms': [s for i, s in enumerate(srs_t)]})
+        # -----------------------------------------------
+        elif acq_type == 'breakdyn':
+            t0_dyn = min(t_starts, key=lambda x: abs(x - 0))
+            t1_dyn = min(t_stops, key=lambda x: abs(x - break_time))
+
+            frm_0 = t_starts.index(t0_dyn)
+            frm_1 = t_stops.index(t1_dyn)
+
+            msrs_class.append({
+                'acq': [acq_type], 'time': (t0_dyn, t1_dyn), 'timings': t_frms,
+                'idxs': (frm_0, frm_1),
+                'frms': [s for i, s in enumerate(srs_t) if i in range(frm_0, frm_1 + 1)]})
+        # -----------------------------------------------
+        elif acq_type == 'fulldyn':
+            t0_dyn = min(t_starts, key=lambda x: abs(x - 0))
+            t1_dyn = min(t_stops, key=lambda x: abs(x - fulldyn_time))
+
+            frm_0 = t_starts.index(t0_dyn)
+            frm_1 = t_stops.index(t1_dyn)
+
+            msrs_class.append({
+                'acq': [acq_type], 'time': (t0_dyn, t1_dyn), 'timings': t_frms,
+                'idxs': (frm_0, frm_1),
+                'frms': [s for i, s in enumerate(srs_t) if i in range(frm_0, frm_1 + 1)]})
+        # -----------------------------------------------
+
+    return {'series': msrs_t, 'descr': msrs_class, 'outpath': amyout}
+
+
+# =====================================================================
+def align_suvr(
+    suvr_tdata,
+    suvr_descr,
+    outpath=None,
+    reg_costfun='nmi',
+    reg_force=False,
+    reg_fwhm=8,
+):
+    '''
+    Align SUVr frames after conversion to NIfTI format.
+
+    Arguments:
+    - reg_constfun: the cost function used in SPM registration/alignment of frames
+    - reg_force:    force running the registration even if the registration results
+                are already calculated and stored in the output folder.
+    - reg_fwhm: the FWHM of the Gaussian kernel used for smoothing the images before
+                registration and only for registration purposes.
+
+    '''
+
+    if outpath is None:
+        align_out = suvr_tdata[next(iter(suvr_tdata))]['files'][0].parent.parent
+    else:
+        align_out = Path(outpath)
+
+    # > NIfTI output folder
+    niidir = align_out / 'NIfTI_SUVr'
+    nimpa.create_dir(niidir)
+
+    # > folder of resampled and aligned NIfTI files (SPM)
+    rsmpl_opth = niidir / 'SPM-aligned'
+    nimpa.create_dir(rsmpl_opth)
+
+    # > the name of the output re-aligned file name
+    faligned = 'SUVr_aligned_' + nimpa.rem_chars(suvr_tdata[next(
+        iter(suvr_tdata))]['series']) + '.nii.gz'
+    faligned = niidir / faligned
+
+    # > Matrices: motion metric + paths to affine
+    R = S = None
+
+    # > check if the file exists
+    if reg_force or not faligned.is_file():
+
+        # > remove any files from previous runs
+        files = niidir.glob('*')
+        for f in files:
+            if f.is_file():
+                os.remove(f)
+            else:
+                shutil.rmtree(f)
+
+        # > output nifty frame files
+        nii_frms = []
+
+        # -----------------------------------------------
+        # > convert the individual DICOM frames to NIfTI
+        for i, k in enumerate(suvr_descr['frms']):
+
+            run([
+                dcm2niix.bin, '-i', 'y', '-v', 'n', '-o', niidir, 'f', '%f_%s',
+                suvr_tdata[k]['files'][0].parent])
+
+            # > get the converted NIfTI file
+            fnii = list(niidir.glob(str(suvr_tdata[k]['tacq']) + '*.nii*'))
+            if len(fnii) != 1:
+                raise ValueError('Unexpected number of converted NIfTI files')
+            else:
+                nii_frms.append(fnii[0])
+        # -----------------------------------------------
+
+        # -----------------------------------------------
+        # > CORE ALIGNMENT OF SUVR FRAMES:
+
+        # > frame-based motion metric (rotations+translation)
+        R = np.zeros((len(nii_frms), len(nii_frms)), dtype=np.float32)
+
+        # > paths to the affine files
+        S = [[None for _ in range(len(nii_frms))] for _ in range(len(nii_frms))]
+
+        # > go through all possible combinations of frame registration
+        for c in combinations(suvr_descr['frms'], 2):
+            frm0 = suvr_descr['frms'].index(c[0])
+            frm1 = suvr_descr['frms'].index(c[1])
+
+            fnii0 = nii_frms[frm0]
+            fnii1 = nii_frms[frm1]
+
+            log.info(f'registration of frame #{frm0} and frame #{frm1}')
+
+            # > one way registration
+            spm_res = nimpa.coreg_spm(fnii0, fnii1, fwhm_ref=reg_fwhm, fwhm_flo=reg_fwhm,
+                                      fwhm=[13, 13], costfun=reg_costfun,
+                                      fcomment=f'_combi_{frm0}-{frm1}', outpath=fnii0.parent,
+                                      visual=0, save_arr=False, del_uncmpr=True)
+
+            S[frm0][frm1] = spm_res['faff']
+
+            rot_ss = np.sum((180 * spm_res['rotations'] / np.pi)**2)**.5
+            trn_ss = np.sum(spm_res['translations']**2)**.5
+            R[frm0, frm1] = rot_ss + trn_ss
+
+            # > the other way registration
+            spm_res = nimpa.coreg_spm(fnii1, fnii0, fwhm_ref=reg_fwhm, fwhm_flo=reg_fwhm,
+                                      fwhm=[13, 13], costfun=reg_costfun,
+                                      fcomment=f'_combi_{frm1}-{frm0}', outpath=fnii0.parent,
+                                      visual=0, save_arr=False, del_uncmpr=True)
+
+            S[frm1][frm0] = spm_res['faff']
+
+            rot_ss = np.sum((180 * spm_res['rotations'] / np.pi)**2)**.5
+            trn_ss = np.sum(spm_res['translations']**2)**.5
+            R[frm1, frm0] = rot_ss + trn_ss
+
+        # > sum frames along floating frames
+        fsum = np.sum(R, axis=0)
+
+        # > sum frames along reference frames
+        rsum = np.sum(R, axis=1)
+
+        # > reference frame for SUVr composite frame
+        rfrm = np.argmin(fsum + rsum)
+
+        niiref = nimpa.getnii(nii_frms[rfrm], output='all')
+
+        # > initialise target aligned SUVr image
+        niiim = np.zeros((len(nii_frms),) + niiref['shape'], dtype=np.float32)
+
+        # > copy in the target frame for SUVr composite
+        niiim[rfrm, ...] = niiref['im']
+
+        for ifrm in range(len(nii_frms)):
+            if ifrm == rfrm:
+                continue
+
+            # > resample images for alignment
+            frsmpl = nimpa.resample_spm(
+                nii_frms[rfrm],
+                nii_frms[ifrm],
+                S[rfrm][ifrm],
+                intrp=1.,
+                outpath=rsmpl_opth,
+                pickname='flo',
+                del_ref_uncmpr=True,
+                del_flo_uncmpr=True,
+                del_out_uncmpr=True,
+            )
+
+            niiim[ifrm, ...] = nimpa.getnii(frsmpl)
+
+        # > save aligned SUVr frames
+        nimpa.array2nii(
+            niiim, niiref['affine'], faligned, descrip='AmyPET: aligned SUVr frames',
+            trnsp=(niiref['transpose'].index(0), niiref['transpose'].index(1),
+                   niiref['transpose'].index(2)), flip=niiref['flip'])
+        # -----------------------------------------------
+
+    return {'fpet': faligned, 'outpath': niidir, 'Metric': R, 'faff': S}
+
+
+# =====================================================================
+
+
+# =====================================================================
+def native_proc(cl_dct, atlas='aal', res='1', outpath=None, refvoi_idx=None, refvoi_name=None):
+    '''
+    Preprocess SPM GM segmentation (from CL output) and AAL atlas
+    to native PET space which is trimmed and upscaled to MR resolution.
+
+    cl_dct:     CL process output dictionary as using the normalisation
+                and segmentation from SPM as used in CL.
+    refvoi_idx: indices of the reference region specific for the chosen
+                atlas.  The reference VOI will then be separately
+                provided in the output.
+    refvoi_name: name of the reference region/VOI
+
+    '''
+
+    # > output path
+    if outpath is None:
+        natout = cl_dct['opth'].parent.parent
+    else:
+        natout = Path(outpath)
+
+    # > get the AAL atlas with the resolution of 1 mm
+    fatl = get_atlas(atlas='aal', res=1, outpath=natout)
+
+    # > trim and upscale the native PET relative to MR resolution
+    trmout = r_trimup(cl_dct['petc']['fim'], cl_dct['mric']['fim'], outpath=natout,
+                      store_img_intrmd=True)
+
+    # > get the trimmed PET as dictionary
+    petdct = nimpa.getnii(trmout['ftrm'], output='all')
+    # > SPM bounding box of the PET image
+    ml = spm12.get_matlab()
+    bbox = spm12.get_bbox(petdct)
+
+    # > get the inverse affine transform to PET native space
+    M = np.linalg.inv(cl_dct['reg2']['affine'])
+    Mm = ml.double(M.tolist())
+
+    # > copy the inverse definitions to be modified with affine to native PET space
+    fmod = shutil.copyfile(cl_dct['norm']['invdef'],
+                           cl_dct['norm']['invdef'].split('.')[0] + '_2nat.nii')
+    eng = spm12.ensure_spm('')
+    eng.amypad_coreg_modify_affine(fmod, Mm)
+
+    # > unzip the atlas and transform it to PET space
+    fniiatl = nimpa.nii_ugzip(fatl, outpath=natout)
+
+    # > inverse transform the atlas to PET space
+    finvatl = spm12.normw_spm(fmod, [fniiatl + ',1'], voxsz=1., intrp=0., bbox=bbox,
+                              outpath=natout)
+
+    # > remove the uncompressed input atlas after transforming it
+    os.remove(fniiatl)
+
+    # > GM mask
+    fgmpet = spm12.resample_spm(trmout['ftrm'], cl_dct['norm']['c1'], M, intrp=1.0, outpath=natout,
+                                pickname='flo', fcomment='_GM_in_PET', del_ref_uncmpr=True,
+                                del_flo_uncmpr=True, del_out_uncmpr=True)
+
+    gm_msk = nimpa.getnii(fgmpet)
+    atl_im = nimpa.getnii(finvatl[0])
+
+    # > get a probability mask for cerebellar GM
+
+    if refvoi_idx is not None:
+        refmsk = np.zeros(petdct['im'].shape, dtype=np.float32)
+        for mi in refvoi_idx:
+            print(mi)
+            msk = atl_im == mi
+            refmsk += msk * gm_msk
+
+    # > probability mask for chosen VOI
+    if refvoi_name is not None:
+        fpmsk = natout / (refvoi_name+'_probmask.nii.gz')
+    else:
+        fpmsk = natout / 'reference_VOI_probmask.nii.gz'
+
+    # > save the mask to NIfTI file
+    nimpa.array2nii(
+        refmsk, petdct['affine'], fpmsk, descrip='AmyPET: probability mask',
+        trnsp=(petdct['transpose'].index(0), petdct['transpose'].index(1),
+               petdct['transpose'].index(2)), flip=petdct['flip'])
+
+    return {
+        'fpet': trmout['ftrm'], 'outpath': natout, 'finvdef': fmod, 'fatl': finvatl, 'fgm': fgmpet,
+        'atlas': atl_im, 'gm_msk': gm_msk, 'frefvoi': fpmsk, 'refvoi': refmsk}
diff --git a/amypet/suvr_tools.py b/amypet/suvr_tools.py
new file mode 100644
index 00000000..22fa13c3
--- /dev/null
+++ b/amypet/suvr_tools.py
@@ -0,0 +1,488 @@
+import logging
+import os
+from pathlib import Path, PurePath
+from subprocess import run
+
+import dcm2niix
+import numpy as np
+from matplotlib import pyplot as plt
+from niftypet import nimpa
+
+logging.basicConfig(level=logging.INFO)
+nifti_ext = ('.nii', '.nii.gz')
+dicom_ext = ('.DCM', '.dcm', '.img', '.IMG', '.ima', '.IMA')
+
+
+# ========================================================================================
+def r_trimup(fpet, fmri, outpath=None, store_img_intrmd=True):
+    '''
+    trim and upscale PET relative to MR T1w or its derivative;
+    derives the scale of upscaling/trimming using the image/voxel sizes
+    '''
+
+    if isinstance(fpet, (str, PurePath)):
+        petdct = nimpa.getnii(fpet, output='all')
+    elif isinstance(fpet, dict) and 'hdr' in fpet:
+        petdct = fpet
+    else:
+        raise ValueError('wrong PET input - accepted are path to image file or dictionary')
+
+    if isinstance(fmri, (str, PurePath)):
+        mridct = nimpa.getnii(fmri, output='all')
+    elif isinstance(fmri, dict) and 'hdr' in fmri:
+        mridct = fmri
+    else:
+        raise ValueError('wrong MR input - accepted are path to image file or dictionary')
+
+    # > get the voxel sizes
+    pet_szyx = petdct['hdr']['pixdim'][1:4]
+    mri_szyx = mridct['hdr']['pixdim'][1:4]
+
+    # > estimate the scale
+    scale = np.abs(np.round(pet_szyx[::-1] / mri_szyx[::-1])).astype(np.int32)
+
+    # > trim the PET image for more accurate regional sampling
+    ftrm = nimpa.imtrimup(fpet, scale=scale, store_img_intrmd=store_img_intrmd, outpath=outpath)
+
+    # > trimmed folder
+    trmdir = Path(ftrm['fimi'][0]).parent
+
+    return {'im': ftrm['im'], 'trmdir': trmdir, 'ftrm': ftrm['fimi'][0], 'trim_scale': scale}
+
+
+# ========================================================================================
+def extract_vois(impet, imlabel, voi_dct, outpath=None, output_masks=False):
+    '''
+    Extract VOI mean values from PET image `impet` using image labels `imlabel`.
+    Both can be dictionaries, file paths or Numpy arrays.
+    They have to be aligned and have the same dimensions.
+    If path (output) is given, the ROI masks will be saved to file(s).
+    Arguments:
+        - impet:    PET image as Numpy array
+        - imlabel:  image of labels (integer values); the labels can come
+                    from T1w-based parcellation or an atlas.
+        - voi_dct:  dictionary of VOIs, with entries of labels creating
+                    composite volumes
+        - output_masks: if `True`, output Numpy VOI masks in the output
+                    dictionary
+        - outpath:  if given as a folder path, the VOI masks will be saved
+    '''
+
+    # > assume none of the below are given
+    # > used only for saving ROI mask to file if requested
+    affine, flip, trnsp = None, None, None
+
+    # ----------------------------------------------
+    # PET
+    if isinstance(impet, dict):
+        im = impet['im']
+        if 'affine' in impet:
+            affine = impet['affine']
+        if 'flip' in impet:
+            flip = impet['flip']
+        if 'transpose' in impet:
+            trnsp = impet['transpose']
+
+    elif isinstance(impet, (str, PurePath)) and os.path.isfile(impet):
+        imd = nimpa.getnii(impet, output='all')
+        im = imd['im']
+        flip = imd['flip']
+        trnsp = imd['transpose']
+
+    elif isinstance(impet, np.ndarray):
+        im = impet
+    # ----------------------------------------------
+
+    # ----------------------------------------------
+    # LABELS
+    if isinstance(imlabel, dict):
+        lbls = imlabel['im']
+        if 'affine' in imlabel and affine is None:
+            affine = imlabel['affine']
+        if 'flip' in imlabel and flip is None:
+            flip = imlabel['flip']
+        if 'transpose' in imlabel and trnsp is None:
+            trnsp = imlabel['transpose']
+
+    elif isinstance(imlabel, (str, PurePath)) and os.path.isfile(imlabel):
+        prd = nimpa.getnii(imlabel, output='all')
+        lbls = prd['im']
+        if affine is None:
+            affine = prd['affine']
+        if flip is None:
+            flip = prd['flip']
+        if trnsp is None:
+            trnsp = prd['transpose']
+
+    elif isinstance(imlabel, np.ndarray):
+        lbls = imlabel
+
+    # > get rid of NaNs if any in the parcellation/label image
+    lbls[np.isnan(lbls)] = 0
+    # ----------------------------------------------
+
+    # ----------------------------------------------
+    # > output dictionary
+    out = {}
+
+    logging.debug('Extracting volumes of interest (VOIs):')
+    for k, voi in enumerate(voi_dct):
+
+        logging.info(f'  VOI: {voi}')
+
+        # > ROI mask
+        rmsk = np.zeros(lbls.shape, dtype=bool)
+        # > number of voxels in the ROI
+        vxsum = 0
+        # > voxel emission sum
+        emsum = 0
+
+        for ri in voi_dct[voi]:
+            logging.debug(f'   label{ri}')
+            rmsk += np.equal(lbls, ri)
+
+        if outpath is not None and not isinstance(imlabel, np.ndarray):
+            nimpa.create_dir(outpath)
+            fvoi = Path(outpath) / (str(voi) + '_mask.nii.gz')
+            nimpa.array2nii(rmsk.astype(np.int8), affine, fvoi,
+                            trnsp=(trnsp.index(0), trnsp.index(1), trnsp.index(2)), flip=flip)
+        else:
+            fvoi = None
+
+        vxsum += np.sum(rmsk)
+        emsum += np.sum(im[rmsk].astype(np.float64))
+
+        out[voi] = {'vox_no': vxsum, 'sum': emsum, 'avg': emsum / vxsum, 'fvoi': fvoi}
+
+        if output_masks:
+            out[voi]['roimsk'] = rmsk
+
+    # ----------------------------------------------
+
+    return out
+
+
+# ========================================================================================
+def preproc_suvr(pet_path, frames=None, outpath=None, fname=None):
+    ''' Prepare the PET image for SUVr analysis.
+        Arguments:
+        - pet_path: path to the folder of DICOM images, or to the NIfTI file
+        - outpath:  output folder path; if not given will assume the parent
+                    folder of the input image
+        - fname:    core name of the static (SUVr) NIfTI file
+        - frames:   list of frames to be used for SUVr processing
+    '''
+
+    if not os.path.exists(pet_path):
+        raise ValueError('The provided path does not exist')
+
+    # > convert the path to Path object
+    pet_path = Path(pet_path)
+
+    # --------------------------------------
+    # > sort out the output folder
+    if outpath is None:
+        petout = pet_path.parent
+    else:
+        petout = Path(outpath)
+
+    nimpa.create_dir(petout)
+
+    if fname is None:
+        fname = nimpa.rem_chars(pet_path.name.split('.')[0]) + '_static.nii.gz'
+    elif not str(fname).endswith(nifti_ext[1]):
+        fname += '.nii.gz'
+    # --------------------------------------
+
+    # > NIfTI case
+    if pet_path.is_file() and str(pet_path).endswith(nifti_ext):
+        logging.info('PET path exists and it is a NIfTI file')
+
+        fpet_nii = pet_path
+
+    # > DICOM case (if any file inside the folder is DICOM)
+    elif pet_path.is_dir() and any([f.suffix in dicom_ext for f in pet_path.glob('*')]):
+
+        # > get the NIfTi images from previous processing
+        fpet_nii = list(petout.glob(pet_path.name + '*.nii*'))
+
+        if not fpet_nii:
+            run([dcm2niix.bin, '-i', 'y', '-v', 'n', '-o', petout, 'f', '%f_%s', pet_path])
+
+        fpet_nii = list(petout.glob(pet_path.name + '*.nii*'))
+
+        # > if cannot find a file it might be due to spaces in folder/file names
+        if not fpet_nii:
+            fpet_nii = list(petout.glob(pet_path.name.replace(' ', '_') + '*.nii*'))
+
+        if not fpet_nii:
+            raise ValueError('No SUVr NIfTI files found')
+        elif len(fpet_nii) > 1:
+            raise ValueError('Too many SUVr NIfTI files found')
+        else:
+            fpet_nii = fpet_nii[0]
+
+    # > read the dynamic image
+    imdct = nimpa.getnii(fpet_nii, output='all')
+
+    # > number of dynamic frames
+    nfrm = imdct['hdr']['dim'][4]
+
+    # > ensure that the frames exist in part of full dynamic image data
+    if frames and nfrm < max(frames):
+        raise ValueError('The selected frames do not exist')
+    elif not frames:
+        frames = np.arange(nfrm)
+
+    logging.info(f'{nfrm} frames have been found in the dynamic image.')
+
+    # ------------------------------------------
+    # > static image file path
+    fstat = petout / fname
+
+    # > check if the static (for SUVr) file already exists
+    if not fstat.is_file():
+
+        if nfrm > 1:
+            imstat = np.sum(imdct['im'][frames, ...], axis=0)
+        else:
+            imstat = np.squeeze(imdct['im'])
+
+        nimpa.array2nii(
+            imstat, imdct['affine'], fstat,
+            trnsp=(imdct['transpose'].index(0), imdct['transpose'].index(1),
+                   imdct['transpose'].index(2)), flip=imdct['flip'])
+
+        logging.info(f'Saved SUVr file image to: {fstat}')
+    # ------------------------------------------
+
+    return {'fpet_nii': fpet_nii, 'fstat': fstat}
+
+
+# ========================================================================================
+# Extract VOI values for SUVr analysis (main function)
+# ========================================================================================
+
+
+def voi_process(petpth, lblpth, t1wpth, voi_dct=None, ref_voi=None, frames=None, fname=None,
+                t1_bias_corr=True, outpath=None, output_masks=True, save_voi_masks=False,
+                qc_plot=True, reg_fwhm_pet=0, reg_fwhm_mri=0, reg_costfun='nmi', reg_fresh=True):
+    ''' Process PET image for VOI extraction using MR-based parcellations.
+        The T1w image and the labels which are based on the image must be
+        in the same image space.
+
+        Arguments:
+        - petpth:   path to the PET NIfTI image
+        - lblpth:   path to the label NIfTI image (parcellations)
+        - t1wpth:   path to the T1w MRI NIfTI image for registration
+        - voi_dct:  dictionary of VOI definitions
+        - ref_voi:  if given and in `voi_dct` it is used as reference region
+                    for calculating SUVr
+        - frames:   select the frames if multi-frame image given;
+                    by default selects all frames
+        - fname:    the core file name for resulting images
+        - t1_bias_corr: it True, performs bias field correction of the T1w image
+        - outpath:  folder path to the output images, including intermediate
+                    images
+        - output_masks: if True, output VOI sampling masks in the output
+                    dictionary
+        - save_voi_masks: if True, saves all the VOI masks to the `masks` folder
+        - qc_plot:  plots the PET images and overlay sampling, and saves it to
+                    a PNG file; requires `output_masks` to be True.
+        - reg_fwhm: FWHMs of the Gaussian filter applied to PET or MRI images
+                    by default 0 mm;
+        - reg_costfun: cost function used in image registration
+        - reg_fresh:runs fresh registration if True, otherwise uses an existing
+                    one if found.
+
+    '''
+
+    # > output dictionary
+    out = {}
+
+    # > make sure the paths are Path objects
+    petpth = Path(petpth)
+    t1wpth = Path(t1wpth)
+    lblpth = Path(lblpth)
+
+    if outpath is None:
+        outpath = petpth.parent
+    else:
+        outpath = Path(outpath)
+
+    out['input'] = {'fpet': petpth, 'ft1w': t1wpth, 'flbl': lblpth}
+
+    if not (petpth.exists() and t1wpth.is_file() and lblpth.is_file()):
+        raise ValueError('One of the three paths to PET, T1w or label image is incorrect.')
+
+    # > if dictionary is not given, the VOI values will be calculated for each unique
+    # > VOI in the label/parcellation image
+    if voi_dct is None:
+        lbl = nimpa.getnii(lblpth)
+        voi_dct = {int(lab): [int(lab)] for lab in np.unique(lbl)}
+
+    if ref_voi is not None and not all([r in voi_dct for r in ref_voi]):
+        raise ValueError('Not all VOIs listed as reference are in the VOI definition dictionary.')
+
+    # > static (SUVr) image preprocessing
+    suvr_preproc = preproc_suvr(petpth, frames=frames,
+                                outpath=outpath / (petpth.name.split('.')[0] + '_suvr'),
+                                fname=fname)
+
+    out.update(suvr_preproc)
+
+    if t1_bias_corr:
+        out['n4'] = nimpa.bias_field_correction(t1wpth, executable='sitk',
+                                                outpath=suvr_preproc['fstat'].parent.parent)
+        fmri = out['n4']['fim']
+    else:
+        fmri = t1wpth
+
+    # --------------------------------------------------
+    # TRIMMING / UPSCALING
+    # > derive the scale of upscaling/trimming using the current
+    # > image/voxel sizes
+    trmout = r_trimup(suvr_preproc['fstat'], lblpth, store_img_intrmd=True)
+
+    # > trimmed folder
+    trmdir = trmout['trmdir']
+
+    # > trimmed and upsampled PET file
+    out['ftrm'] = trmout['ftrm']
+    out['trim_scale'] = trmout['trim_scale']
+    # --------------------------------------------------
+
+    # > - - - - - - - - - - - - - - - - - - - - - - - -
+    # > parcellations in PET space
+    fplbl = trmdir / '{}_Parcellation_in-upsampled-PET.nii.gz'.format(
+        suvr_preproc['fstat'].name.split('.nii')[0])
+
+    if not fplbl.is_file() or reg_fresh:
+
+        logging.info(f'registration with smoothing of {reg_fwhm_pet}, {reg_fwhm_mri} mm'
+                     ' for reference and floating images respectively')
+
+        spm_res = nimpa.coreg_spm(trmout['ftrm'], fmri, fwhm_ref=reg_fwhm_pet,
+                                  fwhm_flo=reg_fwhm_mri, fwhm=[7, 7], costfun=reg_costfun,
+                                  fcomment='', outpath=trmdir, visual=0, save_arr=False,
+                                  del_uncmpr=True)
+
+        flbl_pet = nimpa.resample_spm(
+            trmout['ftrm'],
+            lblpth,
+            spm_res['faff'],
+            outpath=trmdir,
+            intrp=0.,
+            fimout=fplbl,
+            del_ref_uncmpr=True,
+            del_flo_uncmpr=True,
+            del_out_uncmpr=True,
+        )
+
+    out['flbl'] = fplbl
+    # > - - - - - - - - - - - - - - - - - - - - - - - -
+
+    # > get the label image in PET space
+    plbl_dct = nimpa.getnii(fplbl, output='all')
+
+    # > get the sampling output
+    if save_voi_masks:
+        mask_dir = trmdir / 'masks'
+    else:
+        mask_dir = None
+    voival = extract_vois(trmout['im'], plbl_dct, voi_dct, outpath=mask_dir,
+                          output_masks=output_masks)
+
+    # > calculate SUVr if reference regions is given
+    suvrtxt = None
+    if ref_voi is not None:
+
+        suvr = {}
+
+        suvrtxt = ' '
+        for rvoi in ref_voi:
+            ref = voival[rvoi]['avg']
+            suvr[rvoi] = {}
+            for voi in voi_dct:
+                suvr[rvoi][voi] = voival[voi]['avg'] / ref
+
+            # > get the static trimmed image:
+            imsuvr = nimpa.getnii(out['ftrm'], output='all')
+
+            fsuvr = trmdir / 'SUVr_ref-{}_{}'.format(rvoi, suvr_preproc['fstat'].name)
+            # > save SUVr image
+            nimpa.array2nii(
+                imsuvr['im'] / ref, imsuvr['affine'], fsuvr,
+                trnsp=(imsuvr['transpose'].index(0), imsuvr['transpose'].index(1),
+                       imsuvr['transpose'].index(2)), flip=imsuvr['flip'])
+
+            suvr[rvoi]['fsuvr'] = fsuvr
+
+            if 'suvr' in voi_dct:
+                suvrval = suvr[rvoi]['suvr']
+                suvrtxt += f'$SUVR_\\mathrm{{{rvoi}}}=${suvrval:.3f}; '
+
+        out['suvr'] = suvr
+
+    out['vois'] = voival
+
+    # -----------------------------------------
+    # > QC plot
+    if qc_plot and output_masks:
+        showpet = nimpa.imsmooth(trmout['im'].astype(np.float32), voxsize=plbl_dct['voxsize'],
+                                 fwhm=3.)
+
+        def axrange(prf, thrshld, parts):
+            zs = next(x for x, val in enumerate(prf) if val > thrshld)
+            ze = len(prf) - next(x for x, val in enumerate(prf[::-1]) if val > thrshld)
+            # divide the range in parts
+            p = int((ze-zs) / parts)
+            zn = []
+            for k in range(1, parts):
+                zn.append(zs + k*p)
+            return zn
+
+        # z-profile
+        zn = []
+        thrshld = 100
+        zprf = np.sum(voival['neocx']['roimsk'], axis=(1, 2))
+        zn += axrange(zprf, thrshld, 3)
+
+        zprf = np.sum(voival['cblgm']['roimsk'], axis=(1, 2))
+        zn += axrange(zprf, thrshld, 2)
+
+        mskshow = voival['neocx']['roimsk'] + voival['cblgm']['roimsk']
+
+        xn = []
+        xprf = np.sum(mskshow, axis=(0, 1))
+        xn += axrange(xprf, thrshld, 4)
+
+        fig, ax = plt.subplots(2, 3, figsize=(16, 16))
+
+        for ai, zidx in enumerate(zn):
+            msk = mskshow[zidx, ...]
+            impet = showpet[zidx, ...]
+            ax[0][ai].imshow(impet, cmap='magma', vmax=0.9 * impet.max())
+            ax[0][ai].imshow(msk, cmap='gray_r', alpha=0.25)
+            ax[0][ai].xaxis.set_visible(False)
+            ax[0][ai].yaxis.set_visible(False)
+
+        for ai, xidx in enumerate(xn):
+            msk = mskshow[..., xidx]
+            impet = showpet[..., xidx]
+            ax[1][ai].imshow(impet, cmap='magma', vmax=0.9 * impet.max())
+            ax[1][ai].imshow(msk, cmap='gray_r', alpha=0.25)
+            ax[1][ai].xaxis.set_visible(False)
+            ax[1][ai].yaxis.set_visible(False)
+
+        ax[0, 1].text(0, trmout['im'].shape[1] + 10, suvrtxt, fontsize=12)
+
+        plt.tight_layout()
+
+        fqc = trmdir / f'QC_{petpth.name}_Parcellation-over-upsampled-PET.png'
+        plt.savefig(fqc, dpi=300)
+        plt.close('all')
+        out['fqc'] = fqc
+    # -----------------------------------------
+
+    return out
diff --git a/amypet/web.py b/amypet/web.py
index 060340b1..e0b377a8 100644
--- a/amypet/web.py
+++ b/amypet/web.py
@@ -1,4 +1,7 @@
-from streamlit.bootstrap import load_config_options, run
+try:
+    from streamlit.web.bootstrap import load_config_options, run
+except ImportError: # streamlit<1.12
+    from streamlit.bootstrap import load_config_options, run
 
 from . import _backend_web
 
diff --git a/cl_scripts/centi_fbb.py b/cl_scripts/centi_fbb.py
new file mode 100644
index 00000000..8ad3af97
--- /dev/null
+++ b/cl_scripts/centi_fbb.py
@@ -0,0 +1,115 @@
+"""Run centiloid pipeline level-2 for a new tracer FBB"""
+__author__ = "Pawel J Markiewicz"
+__copyright__ = "Copyright 2022"
+
+import pickle
+from pathlib import Path
+from subprocess import run
+
+import dcm2niix
+from niftypet import nimpa
+
+import amypet
+from amypet import centiloid as centi
+
+drv = Path('/data')
+atlases = drv / 'AMYPET' / 'Atlas' / 'CL_2mm'
+opth = drv / 'AMYPET' / 'CL' / 'FBB'
+
+# ----------------------------------------------------------------------
+# > convert to NIfTI as needed
+if False:
+    dirdata = drv / 'AMYPET' / 'CL' / 'FBB'
+    grps = [x for x in dirdata.glob('*') if x.is_dir() if 'E-25' in str(x) or 'YC-10' in str(x)]
+
+    for g in grps:
+
+        print(f'======================\n{g.name}\n======================')
+
+        nout = g / 'NIfTI'
+        nimpa.create_dir(nout)
+        dcmdir = [x for x in g.glob('*') if x.is_dir()]
+
+        for s in dcmdir:
+
+            run([dcm2niix.bin, '-i', 'y', '-v', 'n', '-o', nout, 'f', '%f_%s', s])
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# ELDRL
+dirfbb = drv / 'AMYPET' / 'CL' / 'FBB' / 'E-25_FBB_90110' / 'NIfTI'
+dirpib = drv / 'AMYPET' / 'CL' / 'FBB' / 'E-25_PiB_5070' / 'NIfTI'
+dirmri = drv / 'AMYPET' / 'CL' / 'FBB' / 'E-25_MR' / 'NIfTI'
+ffbbs = sorted(dirfbb.glob('*.nii'))
+fpibs = sorted(dirpib.glob('*.nii'))
+fmris = sorted(dirmri.glob('*.nii'))
+# amypet.im_check_pairs(fpibs, fmris)
+# amypet.im_check_pairs(ffbbs, fmris)
+
+out_pe = centi.run(fpibs, fmris, atlases, tracer='pib', outpath=opth / 'output_pib_e')
+with open(str(opth / 'output_pib_e.pkl'), 'wb') as f:
+    pickle.dump(out_pe, f)
+
+out_fe = centi.run(ffbbs, fmris, atlases, tracer='new', outpath=opth / 'output_fbb_e',
+                   used_saved=True)
+with open(str(opth / 'output_fbb_e.pkl'), 'wb') as f:
+    pickle.dump(out_fe, f)
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# YC
+dirfbb = drv / 'AMYPET' / 'CL' / 'FBB' / 'YC-10_FBB_90110' / 'NIfTI'
+dirpib = drv / 'AMYPET' / 'CL' / 'FBB' / 'YC-10_PiB_5070' / 'NIfTI'
+dirmri = drv / 'AMYPET' / 'CL' / 'FBB' / 'YC-10_MR' / 'NIfTI'
+ffbbs = sorted(dirfbb.glob('*.nii'))
+fpibs = sorted(dirpib.glob('*.nii'))
+fmris = sorted(dirmri.glob('*.nii'))
+# amypet.im_check_pairs(fpibs, fmris)
+# amypet.im_check_pairs(ffbbs, fmris)
+
+out_py = centi.run(fpibs, fmris, atlases, tracer='pib', outpath=opth / 'output_pib_y')
+with open(str(opth / 'output_pib_y.pkl'), 'wb') as f:
+    pickle.dump(out_py, f)
+
+out_fy = centi.run(ffbbs, fmris, atlases, tracer='new', outpath=opth / 'output_fbb_y')
+with open(str(opth / 'output_fbb_y.pkl'), 'wb') as f:
+    pickle.dump(out_fy, f)
+
+# ----------------------------------------------------------------------
+
+# ======================================================================
+# ----------------------------------------------------------------------
+with open(str(opth / 'output_pib_e.pkl'), 'rb') as f:
+    out_pe = pickle.load(f)
+
+with open(str(opth / 'output_fbb_e.pkl'), 'rb') as f:
+    out_fe = pickle.load(f)
+
+with open(str(opth / 'output_pib_y.pkl'), 'rb') as f:
+    out_py = pickle.load(f)
+
+with open(str(opth / 'output_fbb_y.pkl'), 'rb') as f:
+    out_fy = pickle.load(f)
+# ----------------------------------------------------------------------
+
+# combine young and elderly dictionaries
+out_pe.update(out_py)
+out_fe.update(out_fy)
+outp = out_pe
+outf = out_fe
+
+cal = amypet.calib_tracer(outp, outf)
+
+# > save the transformations from SUVr to SUVr_PiB_Calc
+Tsuvr = amypet.save_suvr2pib(cal, 'fbb')
+
+# TEST
+out_t = centi.run(fpibs[0], fmris[0], atlases, tracer='pib', outpath=opth / 'output_test_pib')
+out_tt = centi.run(ffbbs[0], fmris[0], atlases, tracer='fbb', outpath=opth / 'output_test_fbb',
+                   used_saved=True)
+
+with open(opth / 'cal_data.pkl', 'wb') as f:
+    pickle.dump(cal, f)
+
+with open(opth / 'cal_data.pkl', 'rb') as f:
+    cal_fbb = pickle.load(f)
diff --git a/cl_scripts/centi_fbp.py b/cl_scripts/centi_fbp.py
new file mode 100644
index 00000000..7a37a146
--- /dev/null
+++ b/cl_scripts/centi_fbp.py
@@ -0,0 +1,223 @@
+"""Run centiloid pipeline level-2 for a new tracer FBB"""
+__author__ = "Pawel J Markiewicz"
+__copyright__ = "Copyright 2022"
+
+import os
+import pickle
+import re
+import shutil
+import sys
+from pathlib import Path
+from subprocess import run
+
+import dcm2niix
+import numpy as np
+from niftypet import nimpa
+
+import amypet
+from amypet import centiloid as centi
+
+drv = Path('/data')
+atlases = drv / 'AMYPET' / 'Atlas' / 'CL_2mm'
+dirdata = drv / 'AMYPET' / 'CL' / 'FBP'
+opth = dirdata
+
+# ----------------------------------------------------------------------
+# > convert to NIfTI as needed
+if False:
+    intrm_fldr = 'NIfTI_converted_intrm'
+    final_fldr = 'NIfTI'
+
+    # > folder to be corrected (inconsistency detected upsetting the running of AmyPET)
+    dircorr = dirdata / 'Young_Control_PiB' / '411884_PiB3_AC'
+    if (dircorr / '4118844_PiB3_AC').is_dir():
+        os.rename(dircorr / '4118844_PiB3_AC', dircorr / '411884_PiB3_AC')
+
+    grps = [x for x in dirdata.glob('*') if x.is_dir() if 'Elder' in str(x) or 'Young' in str(x)]
+
+    for g in grps:
+        print(f'======================\n{g.name}\n======================')
+
+        dcmall = g.glob('**/*.dcm')
+        dcmdirs = np.unique([d.parent for d in dcmall])
+        # dreps = [d.parent.parent.parent for d in dcmdirs]
+        # # > DICOM folder counter per subject
+        # dcnt = Counter(dreps)§
+
+        # > all DICOMs per subject dictionary
+        sDCM = {}
+        for d in dcmdirs:
+            prnt = d.parent.parent.parent
+            if prnt not in sDCM:
+                sDCM[prnt] = list(d.glob('*.dcm'))
+            else:
+                sDCM[prnt] += list(d.glob('*.dcm'))
+
+        for s in sDCM:
+
+            # ------------------------------------
+            print(f'i> copying files for subject: {s}')
+            tmpdir = opth / 'tmp' / s.name
+            nimpa.create_dir(tmpdir)
+
+            for f in sDCM[s]:
+                shutil.copy(f, tmpdir)
+
+            print('i> convert DICOMs to NIfTI file...')
+            npth = opth / intrm_fldr / g.name / s.name
+            nimpa.create_dir(npth)
+
+            run([dcm2niix.bin, '-v', 'n', '-o', npth, 'f', '%f_%s', tmpdir])
+            # '-i', 'y',
+
+            # > remove the temporary DICOM files
+            shutil.rmtree(tmpdir)
+            # ------------------------------------
+
+            # ------------------------------------
+            # > find all NIfTI files per subject and combine them when needed
+            spth = opth / final_fldr / g.name
+            nimpa.create_dir(spth)
+
+            tmpim = None
+            for fn in npth.glob('*.nii'):
+                ndct = nimpa.getnii(fn, output='all')
+                if tmpim is None:
+                    tmpim = ndct['im']
+                else:
+                    tmpim += ndct['im']
+
+            # > in case of 4D NIfTI image, sum it here
+            if len(tmpim.shape) > 3:
+                tmpim = np.sum(tmpim, axis=0)
+
+            fnm = re.findall(r'\d+[\w\t]*', s.name)[0]
+            fout = spth / (fnm+'.nii.gz')
+
+            nimpa.array2nii(
+                tmpim.astype(np.float32), ndct['affine'], fout,
+                trnsp=(ndct['transpose'].index(0), ndct['transpose'].index(1),
+                       ndct['transpose'].index(2)), flip=ndct['flip'])
+            # ------------------------------------
+
+# # ==============\
+# fnii = Path('/data/AMYPET/CL/FBP/NIfTI_converted/Young_Control_PiB/101437_PiB3_AC/101437_PiB3_AC_unnamed_20120101142541_574880.nii')
+# ndct = nimpa.getnii(fnii, output='all')
+# tmpim = ndct['im']
+# tmpim = np.sum(tmpim, axis=0)
+# fout = fnii.parent/('cc.nii.gz')
+# # ==============\
+
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# ELDRL
+dirfbp = drv / 'AMYPET' / 'CL' / 'FBP' / 'NIfTI' / 'Elder_subject_florbetapir'
+dirpib = drv / 'AMYPET' / 'CL' / 'FBP' / 'NIfTI' / 'Elder_subject_PiB'
+dirmri = drv / 'AMYPET' / 'CL' / 'FBP' / 'NIfTI' / 'Elder_subject_MRI'
+ffbps = sorted(dirfbp.glob('*.nii*'))
+fpibs = sorted(dirpib.glob('*.nii*'))
+fmris = sorted(dirmri.glob('*.nii*'))
+# amypet.im_check_pairs(fpibs, fmris)
+# amypet.im_check_pairs(ffbps, fmris)
+
+out_pe = centi.run(fpibs, fmris, atlases, tracer='pib', outpath=opth / 'output_pib_e')
+with open(str(opth / 'output_pib_e.pkl'), 'wb') as f:
+    pickle.dump(out_pe, f)
+
+out_fe = centi.run(ffbps, fmris, atlases, tracer='new', outpath=opth / 'output_fbp_e',
+                   used_saved=True)
+with open(str(opth / 'output_fbp_e.pkl'), 'wb') as f:
+    pickle.dump(out_fe, f)
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# YC
+dirfbp = drv / 'AMYPET' / 'CL' / 'FBP' / 'NIfTI' / 'Young_Control_florbetapir'
+dirpib = drv / 'AMYPET' / 'CL' / 'FBP' / 'NIfTI' / 'Young_Control_PiB'
+dirmri = drv / 'AMYPET' / 'CL' / 'FBP' / 'NIfTI' / 'Young_Control_MRI'
+ffbps = sorted(dirfbp.glob('*.nii*'))
+fpibs = sorted(dirpib.glob('*.nii*'))
+fmris = sorted(dirmri.glob('*.nii*'))
+# amypet.im_check_pairs(fpibs, fmris)
+# amypet.im_check_pairs(ffbps, fmris)
+
+out_py = centi.run(fpibs, fmris, atlases, tracer='pib', outpath=opth / 'output_pib_y',
+                   used_saved=True)
+with open(str(opth / 'output_pib_y.pkl'), 'wb') as f:
+    pickle.dump(out_py, f)
+
+out_fy = centi.run(ffbps, fmris, atlases, tracer='new', outpath=opth / 'output_fbp_y')
+with open(str(opth / 'output_fbp_y.pkl'), 'wb') as f:
+    pickle.dump(out_fy, f)
+
+# ----------------------------------------------------------------------
+
+sys.exit()
+
+# ======================================================================
+# ----------------------------------------------------------------------
+with open(str(opth / 'output_pib_e.pkl'), 'rb') as f:
+    out_pe = pickle.load(f)
+
+with open(str(opth / 'output_fbp_e.pkl'), 'rb') as f:
+    out_fe = pickle.load(f)
+
+with open(str(opth / 'output_pib_y.pkl'), 'rb') as f:
+    out_py = pickle.load(f)
+
+with open(str(opth / 'output_fbp_y.pkl'), 'rb') as f:
+    out_fy = pickle.load(f)
+# ----------------------------------------------------------------------
+
+# combine young and elderly dictionaries
+out_pe.update(out_py)
+out_fe.update(out_fy)
+outp = out_pe
+outf = out_fe
+
+cal = amypet.calib_tracer(outp, outf)
+
+# > save the transformations from SUVr to SUVr_PiB_Calc
+Tsuvr = amypet.save_suvr2pib(cal, 'fbp')
+
+# TEST
+out_t = centi.run(fpibs[0], fmris[0], atlases, tracer='pib', outpath=opth / 'output_test_pib')
+out_tt = centi.run(ffbbs[0], fmris[0], atlases, tracer='fbb', outpath=opth / 'output_test_fbb',
+                   used_saved=True)
+
+with open(opth / 'cal_data.pkl', 'wb') as f:
+    pickle.dump(cal, f)
+
+with open(opth / 'cal_data.pkl', 'rb') as f:
+    cal_fbp = pickle.load(f)
+
+# TESTS FOR CHECKING CONSISTENCY WITH PROVIDED DATA
+# >>>>>>> 2c5df0134c48bf9fad58528eeac341da8ecc8d6a
+import matplotlib.pyplot as plt
+import openpyxl as xl
+
+info = xl.load_workbook(dirdata / 'Avid_Centiloid_standard_method.xlsx')
+dat = info['Sheet1']
+istart = 2
+
+pid = [i.value for i in dat['A'][istart:]]
+suvr_wc_fbp = np.array([i.value for i in dat['D'][istart:]])
+suvr_wc_pib = [i.value for i in dat['E'][istart:]]
+suvr_wc_cnv = [i.value for i in dat['G'][istart:]]
+cl_wc_fbp = np.array([i.value for i in dat['I'][istart:]])
+cl_wc_pib = [i.value for i in dat['J'][istart:]]
+
+# > index conversion from the xlsx file to the AmyPET indexes (sorted)
+idxs = [pid.index(int(i)) for i in cal['wc']['sbj']]
+suvrf_avid = suvr_wc_fbp[idxs]
+suvrf_amyp = cal['wc']['calib']['cl_suvr'][:, 2]
+
+clf_avid = cl_wc_fbp[idxs]
+
+plt.figure()
+plt.plot(suvrf_avid, suvrf_amyp, '.')
+
+fig, ax = plt.subplots()
+ax.scatter(cal[rvoi]['calib']['cl_suvr'][:, 1], cal[rvoi]['calib']['cl_suvr'][:, 2], c='black')
+amypet.aux.identity_line(ax=ax, ls='--', c='b')
diff --git a/cl_scripts/centi_flt.py b/cl_scripts/centi_flt.py
new file mode 100644
index 00000000..ae5e446d
--- /dev/null
+++ b/cl_scripts/centi_flt.py
@@ -0,0 +1,90 @@
+"""Run centiloid pipeline level-2 for a new tracer FLUTE"""
+__author__ = "Pawel J Markiewicz"
+__copyright__ = "Copyright 2022"
+
+import pickle
+from pathlib import Path
+
+import amypet
+from amypet import centiloid as centi
+
+drv = Path('/data/AMYPET')
+atlases = drv / 'Atlas' / 'CL_2mm'
+opth = drv / 'CL' / 'FLUTE'
+
+# ----------------------------------------------------------------------
+# AD
+dirflt = drv / 'CL' / 'FLUTE' / 'GE_AD_F18_NIFTI'
+dirpib = drv / 'CL' / 'FLUTE' / 'GE_AD_C11_NIFTI'
+dirmri = drv / 'CL' / 'FLUTE' / 'GE_AD_MRI_NIFTI'
+fflts = sorted(dirflt.glob('*.nii'))
+fpibs = sorted(dirpib.glob('*.nii'))
+fmris = sorted(dirmri.glob('*.nii'))
+# amypet.im_check_pairs(fpibs, fmris)
+# amypet.im_check_pairs(fflts, fmris)
+
+out_pe = centi.run(fpibs, fmris, atlases, tracer='pib',
+                   outpath=opth / 'output_pib_a')       # , used_saved=True
+with open(str(opth / 'output_pib_a.pkl'), 'wb') as f:
+    pickle.dump(out_pe, f)
+
+out_fe = centi.run(fflts, fmris, atlases, tracer='new', outpath=opth / 'output_flt_a')
+with open(str(opth / 'output_flt_a.pkl'), 'wb') as f:
+    pickle.dump(out_fe, f)
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# YC
+dirflt = drv / 'CL' / 'FLUTE' / 'GE_YHV_F18_NIFTI'
+dirpib = drv / 'CL' / 'FLUTE' / 'GE_YHV_C11_NIFTI'
+dirmri = drv / 'CL' / 'FLUTE' / 'GE_YHV_MRI_NIFTI'
+fflts = sorted(dirflt.glob('*.nii'))
+fpibs = sorted(dirpib.glob('*.nii'))
+fmris = sorted(dirmri.glob('*.nii'))
+# amypet.im_check_pairs(fpibs, fmris)
+# amypet.im_check_pairs(fflts, fmris)
+
+out_py = centi.run(fpibs, fmris, atlases, tracer='pib', outpath=opth / 'output_pib_y')
+with open(str(opth / 'output_pib_y.pkl'), 'wb') as f:
+    pickle.dump(out_py, f)
+
+out_fy = centi.run(fflts, fmris, atlases, tracer='new', outpath=opth / 'output_flt_y')
+with open(str(opth / 'output_flt_y.pkl'), 'wb') as f:
+    pickle.dump(out_fy, f)
+# ----------------------------------------------------------------------
+
+# ======================================================================
+# ----------------------------------------------------------------------
+with open(str(opth / 'output_pib_a.pkl'), 'rb') as f:
+    out_pa = pickle.load(f)
+
+with open(str(opth / 'output_flt_a.pkl'), 'rb') as f:
+    out_fa = pickle.load(f)
+
+with open(str(opth / 'output_pib_y.pkl'), 'rb') as f:
+    out_py = pickle.load(f)
+
+with open(str(opth / 'output_flt_y.pkl'), 'rb') as f:
+    out_fy = pickle.load(f)
+# ----------------------------------------------------------------------
+
+# combine young and elderly dictionaries
+out_pa.update(out_py)
+out_fa.update(out_fy)
+outp = out_pa
+outf = out_fa
+
+cal = amypet.calib_tracer(outp, outf)
+
+# > save the transformations from SUVr to SUVr_PiB_Calc
+Tsuvr = amypet.save_suvr2pib(cal, 'flute')
+
+# TEST
+out_t = centi.run(fpibs[0], fmris[0], atlases, tracer='pib', outpath=opth / 'output_test_pib')
+out_tt = centi.run(fflts[0], fmris[0], atlases, tracer='flute', outpath=opth / 'output_test_flt')
+
+with open(opth / 'cal_data.pkl', 'wb') as f:
+    pickle.dump(cal, f)
+
+with open(opth / 'cal_data.pkl', 'rb') as f:
+    cal_flt = pickle.load(f)
diff --git a/cl_scripts/centi_pib.py b/cl_scripts/centi_pib.py
new file mode 100644
index 00000000..69219a9b
--- /dev/null
+++ b/cl_scripts/centi_pib.py
@@ -0,0 +1,198 @@
+"""Run calibration centiloid pipeline for PiB"""
+__author__ = "Pawel J Markiewicz"
+__copyright__ = "Copyright 2022"
+import pickle
+from pathlib import Path
+
+import matplotlib.pyplot as plt
+import numpy as np
+import openpyxl as xl
+from scipy.stats import linregress
+
+import amypet
+from amypet.aux import cl_anchor_fldr
+from amypet.centiloid import run as centiloid_run
+
+# > input paths
+drv = Path('/data/AMYPET')
+atlases = drv / 'Atlas' / 'CL_2mm'
+opth = drv / 'CL' / 'PiB'
+
+# ----------------------------------------------------------------------
+# AD
+dirpet = drv / 'CL' / 'PiB' / 'AD-100_PET_5070' / 'nifti'
+dirmri = drv / 'CL' / 'PiB' / 'AD-100_MR' / 'nifti'
+fpets = sorted(dirpet.glob("*.*"))
+fmris = sorted(dirmri.glob("*.*"))
+
+# > run visual check of the images before running the CL pipeline
+# amypet.im_check_pairs(fpets, fmris)
+
+flip_pet = len(fpets) * [(1, 1, 1)]
+flip_pet[9] = (1, -1, 1)
+flip_pet[22] = (1, -1, 1)
+flip_pet[26] = (1, -1, 1)
+flip_pet[34] = (1, -1, 1)
+flip_pet[36] = (1, -1, 1)
+flip_pet[40] = (1, -1, 1)
+flip_pet[44] = (1, -1, 1)
+
+out_ad = centiloid_run(fpets, fmris, atlases, flip_pet=flip_pet, outpath=opth / 'output_pib_ad')
+with open(str(opth / 'output_pib_ad.pkl'), 'wb') as f:
+    pickle.dump(out_ad, f)
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+# YC
+dirpet = drv / 'CL' / 'PiB' / 'YC-0_PET_5070' / 'nifti'
+dirmri = drv / 'CL' / 'PiB' / 'YC-0_MR' / 'nifti'
+fpets = sorted(dirpet.glob("*.*"))
+fmris = sorted(dirmri.glob("*.*"))
+
+# > run visual check of the images before running the CL pipeline
+# amypet.im_check_pairs(fpets, fmris)
+
+out_yc = centiloid_run(fpets, fmris, atlases, outpath=opth / 'output_pib_yc')
+with open(str(opth / 'output_pib_yc.pkl'), 'wb') as f:
+    pickle.dump(out_yc, f)
+# ----------------------------------------------------------------------
+
+# ----------------------------------------------------------------------
+with open(str(opth / 'output_pib_yc.pkl'), 'rb') as f:
+    out_yc = pickle.load(f)
+
+with open(str(opth / 'output_pib_ad.pkl'), 'rb') as f:
+    out_ad = pickle.load(f)
+# ----------------------------------------------------------------------
+
+refs = amypet.get_clref(opth / 'SupplementaryTable1.xlsx')
+
+diff = amypet.check_suvrs(out_yc, out_ad, refs)
+
+diff = amypet.check_cls(out_yc, out_ad, diff, refs)
+
+# cla = amypet.save_cl_anchors(diff, outpath=Path('/home/pawel/Dropbox'))
+cla = amypet.save_cl_anchors(diff)
+
+# > load/test the anchor and conversion tables
+cl_fldr = cl_anchor_fldr
+cla = amypet.get_cl_anchors(cl_fldr)
+cla = amypet.get_cl_anchors(cl_fldr.parent)
+
+tracer = 'flute'
+tracer = 'fbb'
+tracer = 'fbp'
+clc = amypet.get_suvr2pib(tracer, path=cl_fldr)
+clc = amypet.get_suvr2pib(tracer, path=cl_fldr.parent)
+
+
+# ----------------------------------------------------------------------
+def identity_line(ax=None, ls='--', *args, **kwargs):
+    '''
+    plot identity line for any correlation analysis
+    needed for the CL calibration and beyond.
+    '''
+    # see: https://stackoverflow.com/q/22104256/3986320
+    ax = ax or plt.gca()
+    identity, = ax.plot([], [], ls=ls, *args, **kwargs)
+
+    def callback(axes):
+        low_x, high_x = ax.get_xlim()
+        low_y, high_y = ax.get_ylim()
+        low = min(low_x, low_y)
+        high = max(high_x, high_y)
+        identity.set_data([low, high], [low, high])
+
+    callback(ax)
+    ax.callbacks.connect('xlim_changed', callback)
+    ax.callbacks.connect('ylim_changed', callback)
+    return ax
+
+
+# ===================================================
+drv = Path('/Users/pawel/Dropbox/tmp/AMYPAD/AmyPET/')
+
+# - - - - - - - - - - - - - - - - - - - - - - - - -
+with open(drv / 'cal_fbb.pkl', 'rb') as f:
+    cal = pickle.load(f)
+
+info = xl.load_workbook(drv / 'FBBproject_SupplementaryTable.xlsx')
+dat = info['18F-FBB']
+cl_o = np.array([i.value for i in dat['H'][4:29]] + [i.value for i in dat['H'][30:]])
+
+pid = [i.value for i in dat['A'][4:29]] + [int(i.value[1:]) for i in dat['A'][30:]]
+idxs = [pid.index(int(i)) for i in cal['wc']['sbj']]
+cl_os = cl_o[idxs]
+
+rvoi = 'wc'
+cl_amy = cal[rvoi]['calib']['cl_std_fbb']
+
+fig, ax = plt.subplots()
+ax.scatter(cl_os, cl_amy, c='black')
+identity_line(ax=ax, ls='--', c='b')
+ax.set_xlabel('Original F18 CLs')
+ax.set_ylabel('AmyPET F18 CLs')
+ax.grid('on')
+ax.set_title('FBB')
+
+m, a, r, p, stderr = linregress(cl_os, cl_amy)
+r2 = r**2
+ax.text(0, 125, f'$R^2={r2:.4f}$', fontsize=12)
+# - - - - - - - - - - - - - - - - - - - - - - - - -
+
+# - - - - - - - - - - - - - - - - - - - - - - - - -
+with open(drv / 'cal_flt.pkl', 'rb') as f:
+    cal = pickle.load(f)
+
+info = xl.load_workbook(drv / 'flutemetamol.xlsx')
+dat = info['GE PiB & Flutemetamol']
+cl_o = np.array([i.value for i in dat['T'][3:53]] + [i.value for i in dat['T'][54:78]])
+
+pid = [int(i.value[3:]) for i in dat['C'][3:53]] + [int(i.value[4:]) for i in dat['C'][54:78]]
+
+idxs = [pid.index(int(i)) for i in cal['wc']['sbj']]
+cl_os = cl_o[idxs]
+
+rvoi = 'wc'
+cl_amy = cal[rvoi]['calib']['cl_std_fbb']
+
+fig, ax = plt.subplots()
+ax.scatter(cl_os, cl_amy, c='black')
+identity_line(ax=ax, ls='--', c='b')
+ax.set_xlabel('Original F18 CLs')
+ax.set_ylabel('AmyPET F18 CLs')
+ax.grid('on')
+ax.set_title('FLUTE')
+
+m, a, r, p, stderr = linregress(cl_os, cl_amy)
+r2 = r**2
+ax.text(0, 125, f'$R^2={r2:.4f}$', fontsize=12)
+# - - - - - - - - - - - - - - - - - - - - - - - - -
+
+# - - - - - - - - - - - - - - - - - - - - - - - - -
+with open(drv / 'cal_fbp.pkl', 'rb') as f:
+    cal = pickle.load(f)
+
+info = xl.load_workbook(drv / 'Avid_Centiloid_standard_method.xlsx')
+dat = info['Sheet1']
+cl_o = np.array([i.value for i in dat['I'][2:]])
+
+pid = [i.value for i in dat['A'][2:]]
+
+idxs = [pid.index(int(i)) for i in cal['wc']['sbj']]
+cl_os = cl_o[idxs]
+
+rvoi = 'wc'
+cl_amy = cal[rvoi]['calib']['cl_std_fbb']
+
+fig, ax = plt.subplots()
+ax.scatter(cl_os, cl_amy, c='black')
+identity_line(ax=ax, ls='--', c='b')
+ax.set_xlabel('Original F18 CLs')
+ax.set_ylabel('AmyPET F18 CLs')
+ax.grid('on')
+ax.set_title('FBP')
+
+m, a, r, p, stderr = linregress(cl_os, cl_amy)
+r2 = r**2
+ax.text(0, 125, f'$R^2={r2:.4f}$', fontsize=12)
diff --git a/scripts/__init__.py b/scripts/__init__.py
deleted file mode 100644
index e69de29b..00000000
diff --git a/setup.cfg b/setup.cfg
index fc2f522f..6ce5b064 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -29,6 +29,7 @@ classifiers=
     Programming Language :: Python :: 3.7
     Programming Language :: Python :: 3.8
     Programming Language :: Python :: 3.9
+    Programming Language :: Python :: 3.10
     Programming Language :: Python :: 3 :: Only
     Topic :: Scientific/Engineering :: Medical Science Apps.
     Topic :: Software Development :: Libraries
@@ -47,6 +48,7 @@ install_requires=
     setuptools
     shtab>1.3.2
     spm12
+    streamlit>=0.85.2
     tqdm
 include_package_data=True
 packages=find:
@@ -58,7 +60,7 @@ dev=
     pytest-timeout
     pytest-xdist
 gui=Gooey>=1.0.8
-web=streamlit>=0.85.2
+web=
 niftypet=niftypet>=0.0.1
 [options.entry_points]
 console_scripts=
@@ -68,7 +70,7 @@ console_scripts=
 [options.packages.find]
 exclude=tests
 [options.package_data]
-*=*.md, *.rst, *.m, *.png
+*=*.md, *.rst, *.m, *.png, *.nii, *.pkl
 
 [flake8]
 max_line_length=99
@@ -85,6 +87,7 @@ each_dict_entry_on_separate_line=False
 space_between_ending_comma_and_closing_bracket=False
 split_before_named_assigns=False
 split_before_closing_bracket=False
+blank_line_before_nested_class_or_def=0
 
 [isort]
 profile=black


Test patch:
diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index d82d0268..6aa2d1f6 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -2,14 +2,15 @@ name: Test
 on: [push, pull_request]
 jobs:
   check:
-    if: github.event_name != 'pull_request' || github.repository_owner != 'AMYPAD'
+    if: github.event_name != 'pull_request' || !contains('OWNER,MEMBER,COLLABORATOR', github.event.pull_request.author_association)
     runs-on: ubuntu-latest
+    name: Check
     steps:
-    - uses: actions/checkout@v2
-    - uses: actions/setup-python@v2
+    - uses: actions/checkout@v3
+    - uses: actions/setup-python@v4
     - name: set PYSHA
       run: echo "PYSHA=$(python -VV | sha256sum | cut -d' ' -f1)" >> $GITHUB_ENV
-    - uses: actions/cache@v1
+    - uses: actions/cache@v3
       with:
         path: ~/.cache/pre-commit
         key: pre-commit|${{ env.PYSHA }}|${{ hashFiles('.pre-commit-config.yaml') }}
@@ -31,31 +32,31 @@ jobs:
         EVENT: ${{ github.event_name }}
     - run: SKIP=mbeautify pre-commit run -a --show-diff-on-failure
   test:
-    if: github.event_name != 'pull_request' || github.repository_owner != 'AMYPAD'
-    name: py${{ matrix.python }}
+    if: github.event_name != 'pull_request' || !contains('OWNER,MEMBER,COLLABORATOR', github.event.pull_request.author_association)
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python: [3.6, 3.9]
+        python: [3.6, '3.10']
+    name: Test py${{ matrix.python }}
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
       with:
         fetch-depth: 0
-    - uses: actions/setup-python@v2
+    - uses: actions/setup-python@v4
       with:
         python-version: ${{ matrix.python }}
     - run: pip install -U -e .[dev]
     - run: pytest
     - uses: codecov/codecov-action@v1
   cuda:
-    if: github.event_name != 'pull_request' || github.repository_owner != 'AMYPAD'
-    name: CUDA py${{ matrix.python }}
+    if: github.event_name != 'pull_request' || !contains('OWNER,MEMBER,COLLABORATOR', github.event.pull_request.author_association)
     runs-on: [self-hosted, python, matlab, cuda]
     strategy:
       matrix:
         python: [3.9]
+    name: CUDA py${{ matrix.python }}
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
       with:
         fetch-depth: 0
     - name: Run setup-python
@@ -72,10 +73,10 @@ jobs:
     name: PyPI Deploy
     runs-on: ubuntu-latest
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
       with:
         fetch-depth: 0
-    - uses: actions/setup-python@v2
+    - uses: actions/setup-python@v4
     - id: dist
       uses: casperdcl/deploy-pypi@v2
       with: