From d827853b38b941290bdbc9aae1ab2a52d7f88fe4 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sat, 30 Oct 2021 00:25:50 +0100
Subject: [PATCH 01/32] bprj: add output=None convenience

- follow-up to #36
- part of #33
---
 niftypet/nipet/prj/mmrprj.py | 12 ++++++++++--
 1 file changed, 10 insertions(+), 2 deletions(-)

diff --git a/niftypet/nipet/prj/mmrprj.py b/niftypet/nipet/prj/mmrprj.py
index 675edd8f..83e8951f 100644
--- a/niftypet/nipet/prj/mmrprj.py
+++ b/niftypet/nipet/prj/mmrprj.py
@@ -143,7 +143,7 @@ def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=F
 # ------------------------------------------------------------------------
 
 
-def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False):
+def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, output=None):
     '''
     Calculate forward projection for the provided input image.
     Arguments:
@@ -152,6 +152,7 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False):
             transaxial and axial look up tables (LUT).
         isub -- array of transaxial indices of all sinograms (angles x bins) used for subsets;
             when the first element is negative, all transaxial bins are used (as in pure EM-ML).
+        output(CuVec, optional) -- output image.
     '''
     # Get particular scanner parameters: Constants, transaxial and axial LUTs
     Cnt = scanner_params['Cnt']
@@ -199,7 +200,14 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False):
         nvz = Cnt['rSZ_IMZ']
     else:
         nvz = Cnt['SZ_IMZ']
-    bimg = cu.zeros((Cnt['SZ_IMX'], Cnt['SZ_IMY'], nvz), dtype=np.float32)
+
+    out_shape = Cnt['SZ_IMX'], Cnt['SZ_IMY'], nvz
+    if output is None:
+        bimg = cu.zeros(out_shape, dtype=np.float32)
+    else:
+        bimg = cu.asarray(output)
+        assert bimg.shape == out_shape
+        assert bimg.dtype == np.dtype('float32')
 
     # > run back-projection
     petprj.bprj(bimg.cuvec, cu.asarray(sinog).cuvec, txLUT, axLUT, isub, Cnt)

From 1506b963d826b5d3d75dc456787600c1f5488ebc Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sat, 30 Oct 2021 01:16:49 +0100
Subject: [PATCH 02/32] fix tests

- conda py3.6 => 3.7
- matlab py3.9 => 3.8
---
 .github/workflows/test.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 4e93ce37..2593fa48 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -39,7 +39,7 @@ jobs:
     runs-on: [self-hosted, python, cuda, matlab]
     strategy:
       matrix:
-        python: [3.6, 3.9]
+        python: [3.7, 3.8]
     steps:
     - uses: actions/checkout@v2
       with:

From 1da5d5d39fc5457a35e3fffc199c3fa964c97266 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sun, 31 Oct 2021 23:51:20 +0000
Subject: [PATCH 03/32] mmrprj: add `back_prj(div_sino=None)` CUDA convenience

---
 niftypet/nipet/prj/mmrprj.py         | 11 +++++--
 niftypet/nipet/prj/src/prj_module.cu | 19 ++++++++----
 niftypet/nipet/prj/src/prjb.cu       | 45 ++++++++++++++++++----------
 niftypet/nipet/prj/src/prjb.h        |  3 +-
 4 files changed, 54 insertions(+), 24 deletions(-)

diff --git a/niftypet/nipet/prj/mmrprj.py b/niftypet/nipet/prj/mmrprj.py
index 83e8951f..f8eb2eb1 100644
--- a/niftypet/nipet/prj/mmrprj.py
+++ b/niftypet/nipet/prj/mmrprj.py
@@ -143,7 +143,7 @@ def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=F
 # ------------------------------------------------------------------------
 
 
-def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, output=None):
+def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=None, output=None):
     '''
     Calculate forward projection for the provided input image.
     Arguments:
@@ -152,6 +152,7 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, output=None
             transaxial and axial look up tables (LUT).
         isub -- array of transaxial indices of all sinograms (angles x bins) used for subsets;
             when the first element is negative, all transaxial bins are used (as in pure EM-ML).
+        div_sino -- if specificed, backprojects `sino[isub]/div_sino` instead of `sino`.
         output(CuVec, optional) -- output image.
     '''
     # Get particular scanner parameters: Constants, transaxial and axial LUTs
@@ -176,6 +177,10 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, output=None
 
     # > check first the Siemens default sinogram;
     # > for this default shape only full sinograms are expected--no subsets.
+    orig_sino = sino
+    if div_sino is not None:
+        sino = sino[isub, :]
+        div_sino = cu.asarray(div_sino).cuvec
     if len(sino.shape) == 3:
         if sino.shape[0] != nsinos or sino.shape[1] != Cnt['NSANGLES'] or sino.shape[2] != Cnt[
                 'NSBINS']:
@@ -210,7 +215,9 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, output=None
         assert bimg.dtype == np.dtype('float32')
 
     # > run back-projection
-    petprj.bprj(bimg.cuvec, cu.asarray(sinog).cuvec, txLUT, axLUT, isub, Cnt)
+    petprj.bprj(bimg.cuvec,
+                cu.asarray(sinog if div_sino is None else orig_sino).cuvec, txLUT, axLUT, isub,
+                Cnt, div_sino=div_sino)
 
     if not dev_out:
         # > change from GPU optimised image dimensions to the standard Siemens shape
diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index 9309fea7..5530af00 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -29,7 +29,7 @@ Copyrights: 2019
 //--- Available functions
 static PyObject *trnx_prj(PyObject *self, PyObject *args);
 static PyObject *frwd_prj(PyObject *self, PyObject *args);
-static PyObject *back_prj(PyObject *self, PyObject *args);
+static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs);
 static PyObject *osem_rec(PyObject *self, PyObject *args);
 //---
 
@@ -37,7 +37,7 @@ static PyObject *osem_rec(PyObject *self, PyObject *args);
 static PyMethodDef petprj_methods[] = {
     {"tprj", trnx_prj, METH_VARARGS, "Transaxial projector."},
     {"fprj", frwd_prj, METH_VARARGS, "PET forward projector."},
-    {"bprj", back_prj, METH_VARARGS, "PET back projector."},
+    {"bprj", (PyCFunction)back_prj, METH_VARARGS | METH_KEYWORDS, "PET back projector."},
     {"osem", osem_rec, METH_VARARGS, "OSEM reconstruction of PET data."},
     {NULL, NULL, 0, NULL} // Sentinel
 };
@@ -396,7 +396,7 @@ static PyObject *frwd_prj(PyObject *self, PyObject *args) {
 //==============================================================================
 // B A C K   P R O J E C T O R
 //------------------------------------------------------------------------------
-static PyObject *back_prj(PyObject *self, PyObject *args) {
+static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
   // Structure of constants
   Cnst Cnt;
@@ -419,10 +419,16 @@ static PyObject *back_prj(PyObject *self, PyObject *args) {
   // output backprojected image
   PyCuVec<float> *o_bimg;
 
+  // sinogram divisor
+  PyCuVec<float> *o_div_sino = NULL;
+
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   /* Parse the input tuple */
-  if (!PyArg_ParseTuple(args, "OOOOOO", (PyObject **)&o_bimg, (PyObject **)&o_sino, &o_txLUT,
-                        &o_axLUT, &o_subs, &o_mmrcnst))
+
+  static const char *kwds[] = {"bimg", "sino", "txLUT", "axLUT", "subs", "cnst", "div_sino", NULL};
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOOOOO|O", (char **)kwds, (PyObject **)&o_bimg,
+                                   (PyObject **)&o_sino, &o_txLUT, &o_axLUT, &o_subs, &o_mmrcnst,
+                                   (PyObject **)&o_div_sino))
     return NULL;
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -493,6 +499,7 @@ static PyObject *back_prj(PyObject *self, PyObject *args) {
 
     return NULL;
   }
+  if (Py_None == (PyObject *)o_div_sino) o_div_sino = NULL;
 
   int *subs_ = (int *)PyArray_DATA(p_subs);
   short *s2c = (short *)PyArray_DATA(p_s2c);
@@ -535,7 +542,7 @@ static PyObject *back_prj(PyObject *self, PyObject *args) {
 
   //<><><<><><><><><><><><><><><><><><><><><<><><><><<><><><><><><><><><><><><><><><><><<><><><><><><>
   gpu_bprj(o_bimg->vec.data(), o_sino->vec.data(), li2rng, li2sn, li2nos, s2c, aw2ali, crs, subs,
-           Nprj, Naw, N0crs, Cnt);
+           Nprj, Naw, N0crs, Cnt, o_div_sino ? o_div_sino->vec.data() : nullptr);
   //<><><><><><><><><><><>><><><><><><><><><<><><><><<><><><><><><><><><><><><><><><><><<><><><><><><>
 
   // Clean up
diff --git a/niftypet/nipet/prj/src/prjb.cu b/niftypet/nipet/prj/src/prjb.cu
index 8511427b..bf3c0b25 100644
--- a/niftypet/nipet/prj/src/prjb.cu
+++ b/niftypet/nipet/prj/src/prjb.cu
@@ -30,12 +30,18 @@ __global__ void imReduce(float *imr, float *im, int vz0, int nvz) {
 //===============================================================
 
 //**************** DIRECT ***********************************
-__global__ void bprj_drct(const float *sino, float *im, const float *tt, const unsigned char *tv,
-                          const int *subs, const short snno) {
+__global__ void bprj_drct(const float *sino, float *im, float *div_sino, const float *tt,
+                          const unsigned char *tv, const int *subs, const short snno) {
   int ixt = subs[blockIdx.x]; // transaxial indx
   int ixz = threadIdx.x;      // axial (z)
 
-  float bin = sino[c_li2sn[ixz].x + blockIdx.x * snno];
+  float bin;
+  if (div_sino) {
+    // using full sino with subset divisor (divisor may be smaller)
+    bin = sino[c_li2sn[ixz].x + ixt * snno] / div_sino[c_li2sn[ixz].x + blockIdx.x * snno];
+  } else {
+    bin = sino[c_li2sn[ixz].x + blockIdx.x * snno];
+  }
 
   float z = c_li2rng[ixz].x + .5 * SZ_RING;
   int w = (floorf(.5 * SZ_IMZ + SZ_VOXZi * z));
@@ -83,8 +89,9 @@ __global__ void bprj_drct(const float *sino, float *im, const float *tt, const u
 }
 
 //************** OBLIQUE **************************************************
-__global__ void bprj_oblq(const float *sino, float *im, const float *tt, const unsigned char *tv,
-                          const int *subs, const short snno, const int zoff, const short nil2r_c) {
+__global__ void bprj_oblq(const float *sino, float *im, float *div_sino, const float *tt,
+                          const unsigned char *tv, const int *subs, const short snno,
+                          const int zoff, const short nil2r_c) {
 
   int ixz = threadIdx.x + zoff; // axial (z)
 
@@ -92,8 +99,16 @@ __global__ void bprj_oblq(const float *sino, float *im, const float *tt, const u
 
     int ixt = subs[blockIdx.x]; // blockIdx.x is the transaxial bin index
                                 // bin values to be back projected
-    float bin = sino[c_li2sn[ixz].x + snno * blockIdx.x];
-    float bin_ = sino[c_li2sn[ixz].y + snno * blockIdx.x];
+
+    float bin, bin_;
+    if (div_sino) {
+      // using full sino with subset divisor (divisor may be smaller)
+      bin = sino[c_li2sn[ixz].x + snno * ixt] / div_sino[c_li2sn[ixz].x + snno * blockIdx.x];
+      bin_ = sino[c_li2sn[ixz].y + snno * ixt] / div_sino[c_li2sn[ixz].y + snno * blockIdx.x];
+    } else {
+      bin = sino[c_li2sn[ixz].x + snno * blockIdx.x];
+      bin_ = sino[c_li2sn[ixz].y + snno * blockIdx.x];
+    }
 
     //-------------------------------------------------
     /*** accumulation ***/
@@ -188,8 +203,8 @@ __global__ void bprj_oblq(const float *sino, float *im, const float *tt, const u
 
 //--------------------------------------------------------------------------------------------------
 void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt) {
-
+              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
+              float *_d_div_sino) {
   int dev_id;
   cudaGetDevice(&dev_id);
   if (Cnt.LOG <= LOGDEBUG) printf("i> using CUDA device #%d\n", dev_id);
@@ -277,7 +292,7 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
   //-----------------------------------------------------------------------
 
   //============================================================================
-  bprj_drct<<<Nprj, nrng_c>>>(d_sino, d_imf, d_tt, d_tv, d_subs, snno);
+  bprj_drct<<<Nprj, nrng_c>>>(d_sino, d_imf, _d_div_sino, d_tt, d_tv, d_subs, snno);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
@@ -287,11 +302,11 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
   int Nz = ((Noblq + 127) / 128) * 128;
 
   //============================================================================
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, _d_div_sino, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
   HANDLE_ERROR(cudaGetLastError());
 
   zoff += Nz / 2;
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, _d_div_sino, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
@@ -363,19 +378,19 @@ void rec_bprj(float *d_bimg, float *d_sino, int *d_sub, int Nprj, float *d_tt, u
   if (Cnt.LOG <= LOGDEBUG) printf("i> subset back projection (Nprj=%d)... ", Nprj);
 
   //============================================================================
-  bprj_drct<<<Nprj, NRINGS>>>(d_sino, d_bimg, d_tt, d_tv, d_sub, snno);
+  bprj_drct<<<Nprj, NRINGS>>>(d_sino, d_bimg, nullptr, d_tt, d_tv, d_sub, snno);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
   int zoff = NRINGS;
   //============================================================================
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, nullptr, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
   zoff += Nz / 2;
   //============================================================================
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, nullptr, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
diff --git a/niftypet/nipet/prj/src/prjb.h b/niftypet/nipet/prj/src/prjb.h
index d03b4e19..498b79d8 100644
--- a/niftypet/nipet/prj/src/prjb.h
+++ b/niftypet/nipet/prj/src/prjb.h
@@ -8,7 +8,8 @@
 
 // used from Python
 void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt);
+              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
+              float *_d_div_sino = nullptr);
 
 // to be used within CUDA C reconstruction
 void rec_bprj(float *d_bimg, float *d_sino, int *sub, int Nprj,

From 5f0677f5cedf3ceb42fc35a681090a3e5c0f290b Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sun, 31 Oct 2021 23:52:44 +0000
Subject: [PATCH 04/32] mmrprj: add & expose `sync=True` CUDA optimisation

---
 niftypet/nipet/prj/mmrprj.py         | 11 ++++++----
 niftypet/nipet/prj/src/prj_module.cu | 27 ++++++++++++++---------
 niftypet/nipet/prj/src/prjb.cu       | 30 +++++++++++++++-----------
 niftypet/nipet/prj/src/prjb.h        |  2 +-
 niftypet/nipet/prj/src/prjf.cu       | 32 +++++++++++++++++-----------
 niftypet/nipet/prj/src/prjf.h        |  4 ++--
 6 files changed, 64 insertions(+), 42 deletions(-)

diff --git a/niftypet/nipet/prj/mmrprj.py b/niftypet/nipet/prj/mmrprj.py
index f8eb2eb1..b59b1edc 100644
--- a/niftypet/nipet/prj/mmrprj.py
+++ b/niftypet/nipet/prj/mmrprj.py
@@ -44,7 +44,7 @@ def trnx_prj(scanner_params, sino=None, im=None):
 
 
 def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=False,
-             fullsino_out=True, output=None):
+             fullsino_out=True, output=None, sync=True):
     """
     Calculate forward projection (a set of sinograms) for the provided input image.
     Arguments:
@@ -60,6 +60,7 @@ def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=F
             calculations (attenuation=True), the exponential of the negative of the integrated
             mu-values along LOR path is taken at the end.
         output(CuVec, optional) -- output sinogram.
+        sync(bool) -- whether to `cudaDeviceSynchronize()` after.
     """
     # Get particular scanner parameters: Constants, transaxial and axial LUTs
     Cnt = scanner_params['Cnt']
@@ -121,7 +122,7 @@ def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=F
         assert sinog.shape == out_shape
         assert sinog.dtype == np.dtype('float32')
     # --------------------
-    petprj.fprj(sinog.cuvec, cu.asarray(ims).cuvec, txLUT, axLUT, isub, Cnt, att)
+    petprj.fprj(sinog.cuvec, cu.asarray(ims).cuvec, txLUT, axLUT, isub, Cnt, att, sync=sync)
     # --------------------
 
     # get the sinogram bins in a full sinogram if requested
@@ -143,7 +144,8 @@ def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=F
 # ------------------------------------------------------------------------
 
 
-def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=None, output=None):
+def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=None, output=None,
+             sync=True):
     '''
     Calculate forward projection for the provided input image.
     Arguments:
@@ -154,6 +156,7 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=No
             when the first element is negative, all transaxial bins are used (as in pure EM-ML).
         div_sino -- if specificed, backprojects `sino[isub]/div_sino` instead of `sino`.
         output(CuVec, optional) -- output image.
+        sync(bool) -- whether to `cudaDeviceSynchronize()` after.
     '''
     # Get particular scanner parameters: Constants, transaxial and axial LUTs
     Cnt = scanner_params['Cnt']
@@ -217,7 +220,7 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=No
     # > run back-projection
     petprj.bprj(bimg.cuvec,
                 cu.asarray(sinog if div_sino is None else orig_sino).cuvec, txLUT, axLUT, isub,
-                Cnt, div_sino=div_sino)
+                Cnt, div_sino=div_sino, sync=sync)
 
     if not dev_out:
         # > change from GPU optimised image dimensions to the standard Siemens shape
diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index 5530af00..a60a5362 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -28,7 +28,7 @@ Copyrights: 2019
 
 //--- Available functions
 static PyObject *trnx_prj(PyObject *self, PyObject *args);
-static PyObject *frwd_prj(PyObject *self, PyObject *args);
+static PyObject *frwd_prj(PyObject *self, PyObject *args, PyObject *kwargs);
 static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs);
 static PyObject *osem_rec(PyObject *self, PyObject *args);
 //---
@@ -36,7 +36,7 @@ static PyObject *osem_rec(PyObject *self, PyObject *args);
 //> Module Method Table
 static PyMethodDef petprj_methods[] = {
     {"tprj", trnx_prj, METH_VARARGS, "Transaxial projector."},
-    {"fprj", frwd_prj, METH_VARARGS, "PET forward projector."},
+    {"fprj", (PyCFunction)frwd_prj, METH_VARARGS | METH_KEYWORDS, "PET forward projector."},
     {"bprj", (PyCFunction)back_prj, METH_VARARGS | METH_KEYWORDS, "PET back projector."},
     {"osem", osem_rec, METH_VARARGS, "OSEM reconstruction of PET data."},
     {NULL, NULL, 0, NULL} // Sentinel
@@ -229,7 +229,7 @@ static PyObject *trnx_prj(PyObject *self, PyObject *args) {
 // F O R W A R D   P R O J E C T O R
 //------------------------------------------------------------------------------
 
-static PyObject *frwd_prj(PyObject *self, PyObject *args) {
+static PyObject *frwd_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   // Structure of constants
   Cnst Cnt;
 
@@ -254,10 +254,15 @@ static PyObject *frwd_prj(PyObject *self, PyObject *args) {
   // flag for attenuation factors to be found based on mu-map; if 0 normal emission projection is
   // used
   int att;
+
+  bool SYNC = true; // whether to ensure deviceToHost copy on return
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   /* Parse the input tuple */
-  if (!PyArg_ParseTuple(args, "OOOOOOi", (PyObject **)&o_prjout, (PyObject **)&o_im, &o_txLUT,
-                        &o_axLUT, &o_subs, &o_mmrcnst, &att))
+  static const char *kwds[] = {"sino", "im",  "txLUT", "axLUT", "subs",
+                               "cnst", "att", "sync",  NULL};
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOOOOOi|b", (char **)kwds,
+                                   (PyObject **)&o_prjout, (PyObject **)&o_im, &o_txLUT, &o_axLUT,
+                                   &o_subs, &o_mmrcnst, &att, &SYNC))
     return NULL;
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -373,7 +378,7 @@ static PyObject *frwd_prj(PyObject *self, PyObject *args) {
 
   //<><><><><><><<><><><><><><><><><><><><><><><><><<><><><><><><><><><><><><><><><><><><<><><><><><><><><><><>
   gpu_fprj(o_prjout->vec.data(), o_im->vec.data(), li2rng, li2sn, li2nos, s2c, aw2ali, crs, subs,
-           Nprj, Naw, N0crs, Cnt, att);
+           Nprj, Naw, N0crs, Cnt, att, SYNC);
   //<><><><><><><><<><><><><><><><><><><><><><><><><<><><><><><><><><><><><><><><><><><><<><><><><><><><><><><>
 
   // Clean up
@@ -422,13 +427,15 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   // sinogram divisor
   PyCuVec<float> *o_div_sino = NULL;
 
+  bool SYNC = true; // whether to ensure deviceToHost copy on return
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   /* Parse the input tuple */
 
-  static const char *kwds[] = {"bimg", "sino", "txLUT", "axLUT", "subs", "cnst", "div_sino", NULL};
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOOOOO|O", (char **)kwds, (PyObject **)&o_bimg,
+  static const char *kwds[] = {"bimg", "sino",     "txLUT", "axLUT", "subs",
+                               "cnst", "div_sino", "sync",  NULL};
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOOOOO|Ob", (char **)kwds, (PyObject **)&o_bimg,
                                    (PyObject **)&o_sino, &o_txLUT, &o_axLUT, &o_subs, &o_mmrcnst,
-                                   (PyObject **)&o_div_sino))
+                                   (PyObject **)&o_div_sino, &SYNC))
     return NULL;
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -542,7 +549,7 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
   //<><><<><><><><><><><><><><><><><><><><><<><><><><<><><><><><><><><><><><><><><><><><<><><><><><><>
   gpu_bprj(o_bimg->vec.data(), o_sino->vec.data(), li2rng, li2sn, li2nos, s2c, aw2ali, crs, subs,
-           Nprj, Naw, N0crs, Cnt, o_div_sino ? o_div_sino->vec.data() : nullptr);
+           Nprj, Naw, N0crs, Cnt, o_div_sino ? o_div_sino->vec.data() : nullptr, SYNC);
   //<><><><><><><><><><><>><><><><><><><><><<><><><><<><><><><><><><><><><><><><><><><><<><><><><><><>
 
   // Clean up
diff --git a/niftypet/nipet/prj/src/prjb.cu b/niftypet/nipet/prj/src/prjb.cu
index bf3c0b25..5ca9233c 100644
--- a/niftypet/nipet/prj/src/prjb.cu
+++ b/niftypet/nipet/prj/src/prjb.cu
@@ -204,7 +204,7 @@ __global__ void bprj_oblq(const float *sino, float *im, float *div_sino, const f
 //--------------------------------------------------------------------------------------------------
 void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short *s2c,
               int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
-              float *_d_div_sino) {
+              float *_d_div_sino, bool _sync) {
   int dev_id;
   cudaGetDevice(&dev_id);
   if (Cnt.LOG <= LOGDEBUG) printf("i> using CUDA device #%d\n", dev_id);
@@ -281,9 +281,11 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
   cudaMemcpyToSymbol(c_li2nos, li2nos, nil2r_c * sizeof(char));
 
   cudaEvent_t start, stop;
-  cudaEventCreate(&start);
-  cudaEventCreate(&stop);
-  cudaEventRecord(start, 0);
+  if (_sync) {
+    cudaEventCreate(&start);
+    cudaEventCreate(&stop);
+    cudaEventRecord(start, 0);
+  }
 
   if (Cnt.LOG <= LOGDEBUG) printf("i> calculating image through back projection... ");
 
@@ -328,14 +330,18 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
     if (Cnt.LOG <= LOGDEBUG) printf("i> reduced the axial (z) image size to %d\n", nvz);
   }
 
-  cudaEventRecord(stop, 0);
-  cudaEventSynchronize(stop);
-  // cudaDeviceSynchronize();
-  float elapsedTime;
-  cudaEventElapsedTime(&elapsedTime, start, stop);
-  cudaEventDestroy(start);
-  cudaEventDestroy(stop);
-  if (Cnt.LOG <= LOGDEBUG) printf("DONE in %fs.\n", 0.001 * elapsedTime);
+  if (_sync) {
+    cudaEventRecord(stop, 0);
+    cudaEventSynchronize(stop);
+    // cudaDeviceSynchronize();
+    float elapsedTime;
+    cudaEventElapsedTime(&elapsedTime, start, stop);
+    cudaEventDestroy(start);
+    cudaEventDestroy(stop);
+    if (Cnt.LOG <= LOGDEBUG) printf("DONE in %fs.\n", 0.001 * elapsedTime);
+  } else {
+    if (Cnt.LOG <= LOGDEBUG) printf("DONE.\n");
+  }
 
   HANDLE_ERROR(cudaFree(d_tt));
   HANDLE_ERROR(cudaFree(d_tv));
diff --git a/niftypet/nipet/prj/src/prjb.h b/niftypet/nipet/prj/src/prjb.h
index 498b79d8..9618bc40 100644
--- a/niftypet/nipet/prj/src/prjb.h
+++ b/niftypet/nipet/prj/src/prjb.h
@@ -9,7 +9,7 @@
 // used from Python
 void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short *s2c,
               int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
-              float *_d_div_sino = nullptr);
+              float *_d_div_sino = nullptr, bool _sync = true);
 
 // to be used within CUDA C reconstruction
 void rec_bprj(float *d_bimg, float *d_sino, int *sub, int Nprj,
diff --git a/niftypet/nipet/prj/src/prjf.cu b/niftypet/nipet/prj/src/prjf.cu
index 11434bc8..9c937c79 100644
--- a/niftypet/nipet/prj/src/prjf.cu
+++ b/niftypet/nipet/prj/src/prjf.cu
@@ -207,8 +207,8 @@ __global__ void fprj_oblq(float *sino, const float *im, const float *tt, const u
 
 //--------------------------------------------------------------------------------------------------
 void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
-              char att) {
+              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt, char att,
+              bool _sync) {
   int dev_id;
   cudaGetDevice(&dev_id);
   if (Cnt.LOG <= LOGDEBUG) printf("i> using CUDA device #%d\n", dev_id);
@@ -304,9 +304,11 @@ void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2no
   cudaMemcpyToSymbol(c_li2nos, li2nos, nil2r_c * sizeof(char));
 
   cudaEvent_t start, stop;
-  cudaEventCreate(&start);
-  cudaEventCreate(&stop);
-  cudaEventRecord(start, 0);
+  if (_sync) {
+    cudaEventCreate(&start);
+    cudaEventCreate(&stop);
+    cudaEventRecord(start, 0);
+  }
 
   if (Cnt.LOG <= LOGDEBUG) printf("i> calculating sinograms via forward projection...");
 
@@ -333,14 +335,18 @@ void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2no
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
-  cudaEventRecord(stop, 0);
-  cudaEventSynchronize(stop);
-  // cudaDeviceSynchronize();
-  float elapsedTime;
-  cudaEventElapsedTime(&elapsedTime, start, stop);
-  cudaEventDestroy(start);
-  cudaEventDestroy(stop);
-  if (Cnt.LOG <= LOGDEBUG) printf("DONE in %fs.\n", 0.001 * elapsedTime);
+  if (_sync) {
+    cudaEventRecord(stop, 0);
+    cudaEventSynchronize(stop);
+    // cudaDeviceSynchronize();
+    float elapsedTime;
+    cudaEventElapsedTime(&elapsedTime, start, stop);
+    cudaEventDestroy(start);
+    cudaEventDestroy(stop);
+    if (Cnt.LOG <= LOGDEBUG) printf("DONE in %fs.\n", 0.001 * elapsedTime);
+  } else {
+    if (Cnt.LOG <= LOGDEBUG) printf("DONE.\n");
+  }
 
   if (nvz < SZ_IMZ) HANDLE_ERROR(cudaFree(d_im));
   HANDLE_ERROR(cudaFree(d_tt));
diff --git a/niftypet/nipet/prj/src/prjf.h b/niftypet/nipet/prj/src/prjf.h
index a11512cb..38bd525d 100644
--- a/niftypet/nipet/prj/src/prjf.h
+++ b/niftypet/nipet/prj/src/prjf.h
@@ -7,8 +7,8 @@
 #define PRJF_H
 
 void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
-              char att);
+              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt, char att,
+              bool _sync = true);
 
 void rec_fprj(float *d_sino, float *d_img, int *d_sub, int Nprj,
 

From 8fb23a400377a3ff5e2fb212f02a9b1248f9a16c Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Mon, 1 Nov 2021 22:12:52 +0000
Subject: [PATCH 05/32] remove unneeded aw2ali

---
 niftypet/nipet/prj/src/prj_module.cu | 24 +++---------------------
 niftypet/nipet/prj/src/prjb.cu       |  4 ++--
 niftypet/nipet/prj/src/prjb.h        |  4 ++--
 3 files changed, 7 insertions(+), 25 deletions(-)

diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index a60a5362..a800e7e3 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -461,7 +461,6 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   // transaxial sino LUTs:
   PyObject *pd_crs = PyDict_GetItemString(o_txLUT, "crs");
   PyObject *pd_s2c = PyDict_GetItemString(o_txLUT, "s2c");
-  PyObject *pd_aw2ali = PyDict_GetItemString(o_txLUT, "aw2ali");
 
   //-- get the arrays from the dictionaries
   // axLUTs
@@ -474,12 +473,10 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   p_li2rng = (PyArrayObject *)PyArray_FROM_OTF(pd_li2rng, NPY_FLOAT32, NPY_ARRAY_IN_ARRAY);
 
   // sino to crystal, crystals
-  PyArrayObject *p_s2c = NULL, *p_crs = NULL, *p_aw2ali = NULL;
+  PyArrayObject *p_s2c = NULL, *p_crs = NULL;
   p_s2c = (PyArrayObject *)PyArray_FROM_OTF(pd_s2c, NPY_INT16, NPY_ARRAY_IN_ARRAY);
   p_crs = (PyArrayObject *)PyArray_FROM_OTF(pd_crs, NPY_FLOAT32, NPY_ARRAY_IN_ARRAY);
 
-  p_aw2ali = (PyArrayObject *)PyArray_FROM_OTF(pd_aw2ali, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-
   // subsets if using e.g., OSEM
   PyArrayObject *p_subs = NULL;
   p_subs = (PyArrayObject *)PyArray_FROM_OTF(o_subs, NPY_INT32, NPY_ARRAY_IN_ARRAY);
@@ -487,8 +484,7 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
   /* If that didn't work, throw an exception. */
   if (p_li2rno == NULL || p_li2sn == NULL || p_li2sn1 == NULL || p_li2nos == NULL ||
-      p_aw2ali == NULL || p_s2c == NULL || !o_sino || p_crs == NULL || p_subs == NULL ||
-      p_li2rng == NULL || !o_bimg) {
+      p_s2c == NULL || !o_sino || p_crs == NULL || p_subs == NULL || p_li2rng == NULL || !o_bimg) {
     // axLUTs
     Py_XDECREF(p_li2rno);
     Py_XDECREF(p_li2sn);
@@ -496,8 +492,6 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
     Py_XDECREF(p_li2nos);
     Py_XDECREF(p_li2rng);
 
-    // 2D sino LUT
-    Py_XDECREF(p_aw2ali);
     // sino 2 crystals
     Py_XDECREF(p_s2c);
     Py_XDECREF(p_crs);
@@ -510,7 +504,6 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
   int *subs_ = (int *)PyArray_DATA(p_subs);
   short *s2c = (short *)PyArray_DATA(p_s2c);
-  int *aw2ali = (int *)PyArray_DATA(p_aw2ali);
   short *li2sn;
   if (Cnt.SPN == 11) {
     li2sn = (short *)PyArray_DATA(p_li2sn);
@@ -524,7 +517,6 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   int Nprj = PyArray_DIM(p_subs, 0);
   int N0crs = PyArray_DIM(p_crs, 0);
   int N1crs = PyArray_DIM(p_crs, 1);
-  int Naw = PyArray_DIM(p_aw2ali, 0);
 
   int *subs;
   if (subs_[0] == -1) {
@@ -558,7 +550,6 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   Py_DECREF(p_li2sn);
   Py_DECREF(p_li2sn1);
   Py_DECREF(p_li2nos);
-  Py_DECREF(p_aw2ali);
   Py_DECREF(p_s2c);
   Py_DECREF(p_crs);
   Py_DECREF(p_subs);
@@ -633,7 +624,6 @@ static PyObject *osem_rec(PyObject *self, PyObject *args) {
   // transaxial sino LUTs:
   PyObject *pd_crs = PyDict_GetItemString(o_txLUT, "crs");
   PyObject *pd_s2c = PyDict_GetItemString(o_txLUT, "s2c");
-  PyObject *pd_aw2ali = PyDict_GetItemString(o_txLUT, "aw2ali");
 
   //-- get the arrays from the dictionaries
   // output back-projection image
@@ -673,10 +663,6 @@ static PyObject *osem_rec(PyObject *self, PyObject *args) {
   p_li2nos = (PyArrayObject *)PyArray_FROM_OTF(pd_li2nos, NPY_INT8, NPY_ARRAY_IN_ARRAY);
   p_li2rng = (PyArrayObject *)PyArray_FROM_OTF(pd_li2rng, NPY_FLOAT32, NPY_ARRAY_IN_ARRAY);
 
-  // 2D sino index LUT:
-  PyArrayObject *p_aw2ali = NULL;
-  p_aw2ali = (PyArrayObject *)PyArray_FROM_OTF(pd_aw2ali, NPY_INT32, NPY_ARRAY_IN_ARRAY);
-
   // sino to crystal, crystals
   PyArrayObject *p_s2c = NULL, *p_crs = NULL;
   p_s2c = (PyArrayObject *)PyArray_FROM_OTF(pd_s2c, NPY_INT16, NPY_ARRAY_IN_ARRAY);
@@ -687,7 +673,7 @@ static PyObject *osem_rec(PyObject *self, PyObject *args) {
   if (p_imgout == NULL || p_rcnmsk == NULL || p_subs == NULL || p_psng == NULL || p_rsng == NULL ||
       p_ssng == NULL || p_nsng == NULL || p_asng == NULL || p_imgsens == NULL ||
       p_li2rno == NULL || p_li2sn == NULL || p_li2sn1 == NULL || p_li2nos == NULL ||
-      p_aw2ali == NULL || p_s2c == NULL || p_crs == NULL || p_krnl == NULL) {
+      p_s2c == NULL || p_crs == NULL || p_krnl == NULL) {
     //> output image
     PyArray_DiscardWritebackIfCopy(p_imgout);
     Py_XDECREF(p_imgout);
@@ -713,8 +699,6 @@ static PyObject *osem_rec(PyObject *self, PyObject *args) {
     Py_XDECREF(p_li2sn);
     Py_XDECREF(p_li2sn1);
     Py_XDECREF(p_li2nos);
-    //> 2D sinogram LUT
-    Py_XDECREF(p_aw2ali);
     //> sinogram to crystal LUTs
     Py_XDECREF(p_s2c);
     Py_XDECREF(p_crs);
@@ -757,7 +741,6 @@ static PyObject *osem_rec(PyObject *self, PyObject *args) {
   float *li2rng = (float *)PyArray_DATA(p_li2rng);
   float *crs = (float *)PyArray_DATA(p_crs);
   short *s2c = (short *)PyArray_DATA(p_s2c);
-  int *aw2ali = (int *)PyArray_DATA(p_aw2ali);
 
   int N0crs = PyArray_DIM(p_crs, 0);
   int N1crs = PyArray_DIM(p_crs, 1);
@@ -801,7 +784,6 @@ static PyObject *osem_rec(PyObject *self, PyObject *args) {
   Py_DECREF(p_li2sn);
   Py_DECREF(p_li2sn1);
   Py_DECREF(p_li2nos);
-  Py_DECREF(p_aw2ali);
   Py_DECREF(p_s2c);
   Py_DECREF(p_crs);
 
diff --git a/niftypet/nipet/prj/src/prjb.cu b/niftypet/nipet/prj/src/prjb.cu
index 5ca9233c..5d023a5a 100644
--- a/niftypet/nipet/prj/src/prjb.cu
+++ b/niftypet/nipet/prj/src/prjb.cu
@@ -203,8 +203,8 @@ __global__ void bprj_oblq(const float *sino, float *im, float *div_sino, const f
 
 //--------------------------------------------------------------------------------------------------
 void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
-              float *_d_div_sino, bool _sync) {
+              float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt, float *_d_div_sino,
+              bool _sync) {
   int dev_id;
   cudaGetDevice(&dev_id);
   if (Cnt.LOG <= LOGDEBUG) printf("i> using CUDA device #%d\n", dev_id);
diff --git a/niftypet/nipet/prj/src/prjb.h b/niftypet/nipet/prj/src/prjb.h
index 9618bc40..447392bb 100644
--- a/niftypet/nipet/prj/src/prjb.h
+++ b/niftypet/nipet/prj/src/prjb.h
@@ -7,8 +7,8 @@
 #define PRJB_H
 
 // used from Python
-void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt,
+void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short2 *d_s2c,
+              float4 *d_crs, int *d_subs, float *d_tt, unsigned char *d_tv, int Nprj, Cnst Cnt,
               float *_d_div_sino = nullptr, bool _sync = true);
 
 // to be used within CUDA C reconstruction

From c9c436934c2032086c1b1263620818ef295d342d Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Mon, 1 Nov 2021 22:32:33 +0000
Subject: [PATCH 06/32] clean .cuvec usages

- thanks to https://github.com/AMYPAD/CuVec/pull/17
---
 niftypet/nipet/prj/mmrprj.py         |  7 +++----
 niftypet/nipet/prj/mmrrec.py         |  6 ++----
 niftypet/nipet/prj/mmrsim.py         |  3 +--
 niftypet/nipet/prj/src/prj_module.cu | 21 ++++++++++-----------
 niftypet/nipet/sct/mmrsct.py         |  3 +--
 pyproject.toml                       |  2 +-
 setup.cfg                            |  4 ++--
 7 files changed, 20 insertions(+), 26 deletions(-)

diff --git a/niftypet/nipet/prj/mmrprj.py b/niftypet/nipet/prj/mmrprj.py
index b59b1edc..2244d074 100644
--- a/niftypet/nipet/prj/mmrprj.py
+++ b/niftypet/nipet/prj/mmrprj.py
@@ -122,7 +122,7 @@ def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=F
         assert sinog.shape == out_shape
         assert sinog.dtype == np.dtype('float32')
     # --------------------
-    petprj.fprj(sinog.cuvec, cu.asarray(ims).cuvec, txLUT, axLUT, isub, Cnt, att, sync=sync)
+    petprj.fprj(sinog, cu.asarray(ims), txLUT, axLUT, isub, Cnt, att, sync=sync)
     # --------------------
 
     # get the sinogram bins in a full sinogram if requested
@@ -183,7 +183,7 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=No
     orig_sino = sino
     if div_sino is not None:
         sino = sino[isub, :]
-        div_sino = cu.asarray(div_sino).cuvec
+        div_sino = cu.asarray(div_sino)
     if len(sino.shape) == 3:
         if sino.shape[0] != nsinos or sino.shape[1] != Cnt['NSANGLES'] or sino.shape[2] != Cnt[
                 'NSBINS']:
@@ -218,8 +218,7 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=No
         assert bimg.dtype == np.dtype('float32')
 
     # > run back-projection
-    petprj.bprj(bimg.cuvec,
-                cu.asarray(sinog if div_sino is None else orig_sino).cuvec, txLUT, axLUT, isub,
+    petprj.bprj(bimg, cu.asarray(sinog if div_sino is None else orig_sino), txLUT, axLUT, isub,
                 Cnt, div_sino=div_sino, sync=sync)
 
     if not dev_out:
diff --git a/niftypet/nipet/prj/mmrrec.py b/niftypet/nipet/prj/mmrrec.py
index 0f9d0d57..a2b2411d 100644
--- a/niftypet/nipet/prj/mmrrec.py
+++ b/niftypet/nipet/prj/mmrrec.py
@@ -232,8 +232,7 @@ def osemone(datain, mumaps, hst, scanner_params, recmod=3, itr=4, fwhm=0., psf=N
             log.info('using provided attenuation factor sinogram')
         else:
             asng = cu.zeros(psng.shape, dtype=np.float32)
-            petprj.fprj(asng.cuvec,
-                        cu.asarray(mus).cuvec, txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt,
+            petprj.fprj(asng, cu.asarray(mus), txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt,
                         1)
     # > combine attenuation and normalisation
     ansng = asng * nsng
@@ -294,8 +293,7 @@ def osemone(datain, mumaps, hst, scanner_params, recmod=3, itr=4, fwhm=0., psf=N
         sinoTIdx[n, 0] = Nprj
         sinoTIdx[n, 1:], s = get_subsets14(n, scanner_params)
         # sensitivity image
-        petprj.bprj(tmpsens.cuvec,
-                    cu.asarray(ansng[sinoTIdx[n, 1:], :]).cuvec, txLUT, axLUT, sinoTIdx[n, 1:],
+        petprj.bprj(tmpsens, cu.asarray(ansng[sinoTIdx[n, 1:], :]), txLUT, axLUT, sinoTIdx[n, 1:],
                     Cnt)
         imgsens[n] = tmpsens
     del tmpsens
diff --git a/niftypet/nipet/prj/mmrsim.py b/niftypet/nipet/prj/mmrsim.py
index 2332fa85..1ae5eb3f 100644
--- a/niftypet/nipet/prj/mmrsim.py
+++ b/niftypet/nipet/prj/mmrsim.py
@@ -259,8 +259,7 @@ def simulate_recon(
             sinoTIdx[n, 1:], s = mmrrec.get_subsets14(n, scanner_params)
 
             # > sensitivity image
-            petprj.bprj(tmpsim.cuvec,
-                        cu.asarray(attsino[sinoTIdx[n, 1:], :]).cuvec, txLUT, axLUT,
+            petprj.bprj(tmpsim, cu.asarray(attsino[sinoTIdx[n, 1:], :]), txLUT, axLUT,
                         sinoTIdx[n, 1:], Cnt)
             sim[n] = tmpsim
         del tmpsim
diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index a800e7e3..5f12ccb1 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -243,13 +243,13 @@ static PyObject *frwd_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   PyObject *o_txLUT;
 
   // input image to be forward projected  (reshaped for GPU execution)
-  PyCuVec<float> *o_im;
+  PyCuVec<float> *o_im = NULL;
 
   // subsets for OSEM, first the default
   PyObject *o_subs;
 
   // output projection sino
-  PyCuVec<float> *o_prjout;
+  PyCuVec<float> *o_prjout = NULL;
 
   // flag for attenuation factors to be found based on mu-map; if 0 normal emission projection is
   // used
@@ -260,9 +260,9 @@ static PyObject *frwd_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   /* Parse the input tuple */
   static const char *kwds[] = {"sino", "im",  "txLUT", "axLUT", "subs",
                                "cnst", "att", "sync",  NULL};
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOOOOOi|b", (char **)kwds,
-                                   (PyObject **)&o_prjout, (PyObject **)&o_im, &o_txLUT, &o_axLUT,
-                                   &o_subs, &o_mmrcnst, &att, &SYNC))
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&O&OOOOi|b", (char **)kwds, &asPyCuVec_f,
+                                   &o_prjout, &asPyCuVec_f, &o_im, &o_txLUT, &o_axLUT, &o_subs,
+                                   &o_mmrcnst, &att, &SYNC))
     return NULL;
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -416,13 +416,13 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   PyObject *o_txLUT;
 
   // sino to be back projected to image (both reshaped for GPU execution)
-  PyCuVec<float> *o_sino;
+  PyCuVec<float> *o_sino = NULL;
 
   // subsets for OSEM, first the default
   PyObject *o_subs;
 
   // output backprojected image
-  PyCuVec<float> *o_bimg;
+  PyCuVec<float> *o_bimg = NULL;
 
   // sinogram divisor
   PyCuVec<float> *o_div_sino = NULL;
@@ -433,9 +433,9 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
   static const char *kwds[] = {"bimg", "sino",     "txLUT", "axLUT", "subs",
                                "cnst", "div_sino", "sync",  NULL};
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "OOOOOO|Ob", (char **)kwds, (PyObject **)&o_bimg,
-                                   (PyObject **)&o_sino, &o_txLUT, &o_axLUT, &o_subs, &o_mmrcnst,
-                                   (PyObject **)&o_div_sino, &SYNC))
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&O&OOOO|O&b", (char **)kwds, &asPyCuVec_f,
+                                   &o_bimg, &asPyCuVec_f, &o_sino, &o_txLUT, &o_axLUT, &o_subs,
+                                   &o_mmrcnst, &asPyCuVec_f, &o_div_sino, &SYNC))
     return NULL;
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -500,7 +500,6 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
     return NULL;
   }
-  if (Py_None == (PyObject *)o_div_sino) o_div_sino = NULL;
 
   int *subs_ = (int *)PyArray_DATA(p_subs);
   short *s2c = (short *)PyArray_DATA(p_s2c);
diff --git a/niftypet/nipet/sct/mmrsct.py b/niftypet/nipet/sct/mmrsct.py
index 3245b3bc..de830b48 100644
--- a/niftypet/nipet/sct/mmrsct.py
+++ b/niftypet/nipet/sct/mmrsct.py
@@ -571,8 +571,7 @@ def vsm(
     currentspan = Cnt['SPN']
     Cnt['SPN'] = 1
     atto = cu.zeros((txLUT['Naw'], Cnt['NSN1']), dtype=np.float32)
-    petprj.fprj(atto.cuvec,
-                cu.asarray(mu_sctonly).cuvec, txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt, 1)
+    petprj.fprj(atto, cu.asarray(mu_sctonly), txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt, 1)
     atto = mmraux.putgaps(atto, txLUT, Cnt)
     # --------------------------------------------------------------
     # > get norm components setting the geometry and axial to ones
diff --git a/pyproject.toml b/pyproject.toml
index 786e72f6..93efe209 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [build-system]
 requires = ["setuptools>=42", "wheel", "setuptools_scm[toml]>=3.4",
-            "cuvec>=2.5.0", "ninst>=0.10.0", "numpy>=1.14", "miutil[cuda]>=0.4.0",
+            "cuvec>=2.8.0", "ninst>=0.10.0", "numpy>=1.14", "miutil[cuda]>=0.4.0",
             "scikit-build>=0.11.0", "cmake>=3.18", "ninja"]
 
 [tool.setuptools_scm]
diff --git a/setup.cfg b/setup.cfg
index 06f66b3e..c538bb46 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -39,7 +39,7 @@ setup_requires=
     setuptools>=42
     wheel
     setuptools_scm[toml]
-    cuvec>=2.5.0
+    cuvec>=2.8.0
     miutil[cuda]>=0.4.0
     ninst>=0.10.0
     numpy>=1.14
@@ -47,7 +47,7 @@ setup_requires=
     cmake>=3.18
     ninja
 install_requires=
-    cuvec>=2.5.0
+    cuvec>=2.8.0
     miutil>=0.6.0
     nibabel>=2.4.0
     nimpa>=2.0.0

From 22f8999cc651af6e650e8e30a3523a300ee4f5bf Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Mon, 1 Nov 2021 23:29:00 +0000
Subject: [PATCH 07/32] bubble up some mallocs

---
 niftypet/nipet/prj/src/prj_module.cu | 31 +++++++++++++++++++++++--
 niftypet/nipet/prj/src/prjb.cu       | 34 +++-------------------------
 2 files changed, 32 insertions(+), 33 deletions(-)

diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index 5f12ccb1..ca05a0c8 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -539,8 +539,35 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   HANDLE_ERROR(cudaSetDevice(Cnt.DEVID));
 
   //<><><<><><><><><><><><><><><><><><><><><<><><><><<><><><><><><><><><><><><><><><><><<><><><><><><>
-  gpu_bprj(o_bimg->vec.data(), o_sino->vec.data(), li2rng, li2sn, li2nos, s2c, aw2ali, crs, subs,
-           Nprj, Naw, N0crs, Cnt, o_div_sino ? o_div_sino->vec.data() : nullptr, SYNC);
+  float4 *d_crs;
+  HANDLE_ERROR(cudaMalloc(&d_crs, N0crs * sizeof(float4)));
+  HANDLE_ERROR(cudaMemcpy(d_crs, crs, N0crs * sizeof(float4), cudaMemcpyHostToDevice));
+
+  short2 *d_s2c;
+  HANDLE_ERROR(cudaMalloc(&d_s2c, AW * sizeof(short2)));
+  HANDLE_ERROR(cudaMemcpy(d_s2c, s2c, AW * sizeof(short2), cudaMemcpyHostToDevice));
+
+  float *d_tt;
+  HANDLE_ERROR(cudaMalloc(&d_tt, N_TT * AW * sizeof(float)));
+
+  unsigned char *d_tv;
+  HANDLE_ERROR(cudaMalloc(&d_tv, N_TV * AW * sizeof(unsigned char)));
+  HANDLE_ERROR(cudaMemset(d_tv, 0, N_TV * AW * sizeof(unsigned char)));
+
+  // array of subset projection bins
+  int *d_subs;
+  HANDLE_ERROR(cudaMalloc(&d_subs, Nprj * sizeof(int)));
+  HANDLE_ERROR(cudaMemcpy(d_subs, subs, Nprj * sizeof(int), cudaMemcpyHostToDevice));
+
+  gpu_bprj(o_bimg->vec.data(), o_sino->vec.data(), li2rng, li2sn, li2nos, (short2 *)d_s2c,
+           (float4 *)d_crs, d_subs, d_tt, d_tv, Nprj, Cnt,
+           o_div_sino ? o_div_sino->vec.data() : nullptr, SYNC);
+
+  HANDLE_ERROR(cudaFree(d_subs));
+  HANDLE_ERROR(cudaFree(d_tv));
+  HANDLE_ERROR(cudaFree(d_tt));
+  HANDLE_ERROR(cudaFree(d_s2c));
+  HANDLE_ERROR(cudaFree(d_crs));
   //<><><><><><><><><><><>><><><><><><><><><<><><><><<><><><><><><><><><><><><><><><><><<><><><><><><>
 
   // Clean up
diff --git a/niftypet/nipet/prj/src/prjb.cu b/niftypet/nipet/prj/src/prjb.cu
index 5d023a5a..2af74ca8 100644
--- a/niftypet/nipet/prj/src/prjb.cu
+++ b/niftypet/nipet/prj/src/prjb.cu
@@ -202,35 +202,13 @@ __global__ void bprj_oblq(const float *sino, float *im, float *div_sino, const f
 }
 
 //--------------------------------------------------------------------------------------------------
-void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt, float *_d_div_sino,
-              bool _sync) {
+void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short2 *d_s2c,
+              float4 *d_crs, int *d_subs, float *d_tt, unsigned char *d_tv, int Nprj, Cnst Cnt,
+              float *_d_div_sino, bool _sync) {
   int dev_id;
   cudaGetDevice(&dev_id);
   if (Cnt.LOG <= LOGDEBUG) printf("i> using CUDA device #%d\n", dev_id);
 
-  //--- TRANSAXIAL COMPONENT
-  float4 *d_crs;
-  HANDLE_ERROR(cudaMalloc(&d_crs, N0crs * sizeof(float4)));
-  HANDLE_ERROR(cudaMemcpy(d_crs, crs, N0crs * sizeof(float4), cudaMemcpyHostToDevice));
-
-  short2 *d_s2c;
-  HANDLE_ERROR(cudaMalloc(&d_s2c, AW * sizeof(short2)));
-  HANDLE_ERROR(cudaMemcpy(d_s2c, s2c, AW * sizeof(short2), cudaMemcpyHostToDevice));
-
-  float *d_tt;
-  HANDLE_ERROR(cudaMalloc(&d_tt, N_TT * AW * sizeof(float)));
-
-  unsigned char *d_tv;
-  HANDLE_ERROR(cudaMalloc(&d_tv, N_TV * AW * sizeof(unsigned char)));
-  HANDLE_ERROR(cudaMemset(d_tv, 0, N_TV * AW * sizeof(unsigned char)));
-
-  // array of subset projection bins
-  int *d_subs;
-  HANDLE_ERROR(cudaMalloc(&d_subs, Nprj * sizeof(int)));
-  HANDLE_ERROR(cudaMemcpy(d_subs, subs, Nprj * sizeof(int), cudaMemcpyHostToDevice));
-  //---
-
   //-----------------------------------------------------------------
   // RINGS: either all or a subset of rings can be used for fast calc.
   //-----------------------------------------------------------------
@@ -342,12 +320,6 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
   } else {
     if (Cnt.LOG <= LOGDEBUG) printf("DONE.\n");
   }
-
-  HANDLE_ERROR(cudaFree(d_tt));
-  HANDLE_ERROR(cudaFree(d_tv));
-  HANDLE_ERROR(cudaFree(d_subs));
-  HANDLE_ERROR(cudaFree(d_crs));
-  HANDLE_ERROR(cudaFree(d_s2c));
 }
 
 //=======================================================================

From dc69d83d3f112b0cd41422716092e30b2bbaf586 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Tue, 2 Nov 2021 00:02:52 +0000
Subject: [PATCH 08/32] mmraux.remgaps: use CuVec

---
 niftypet/nipet/CMakeLists.txt    |  2 +
 niftypet/nipet/mmraux.py         | 16 +++++---
 niftypet/nipet/src/aux_module.cu | 66 ++++++++------------------------
 3 files changed, 29 insertions(+), 55 deletions(-)

diff --git a/niftypet/nipet/CMakeLists.txt b/niftypet/nipet/CMakeLists.txt
index 3dac1042..70b59188 100644
--- a/niftypet/nipet/CMakeLists.txt
+++ b/niftypet/nipet/CMakeLists.txt
@@ -6,6 +6,7 @@ project(mmr_auxe)
 file(GLOB SRC LIST_DIRECTORIES false "src/*.cu")
 include_directories(src)
 include_directories(${Python3_INCLUDE_DIRS})
+include_directories(${CUVEC_INCLUDE_DIRS})
 include_directories(${Python3_NumPy_INCLUDE_DIRS})
 
 add_library(${PROJECT_NAME} SHARED ${SRC})
@@ -19,6 +20,7 @@ if(SKBUILD)
 python_extension_module(${PROJECT_NAME})
 endif()
 set_target_properties(${PROJECT_NAME} PROPERTIES
+  CXX_STANDARD 11
   #VERSION ${CMAKE_PROJECT_VERSION} SOVERSION ${CMAKE_PROJECT_VERSION_MAJOR}
   PREFIX ""  # remove shared lib prefix to make importable
   INTERFACE_${PROJECT_NAME}_MAJOR_VERSION ${CMAKE_PROJECT_VERSION_MAJOR})
diff --git a/niftypet/nipet/mmraux.py b/niftypet/nipet/mmraux.py
index 60b45d0d..a70f4f13 100644
--- a/niftypet/nipet/mmraux.py
+++ b/niftypet/nipet/mmraux.py
@@ -10,6 +10,7 @@
 from pathlib import Path
 from textwrap import dedent
 
+import cuvec as cu
 import numpy as np
 import pydicom as dcm
 from miutil.fdio import hasext
@@ -1142,18 +1143,21 @@ def putgaps(s, txLUT, Cnt, sino_no=0):
     return sino.astype(s.dtype)
 
 
-def remgaps(sino, txLUT, Cnt):
-
+def remgaps(sino, txLUT, Cnt, output=None):
     # number of sino planes (2D sinos) depends on the span used
     nsinos = sino.shape[0]
 
-    # preallocate output sino without gaps, always in float
-    s = np.zeros((txLUT['Naw'], nsinos), dtype=np.float32)
+    # preallocate output sino without gaps
+    if output is None:
+        output = cu.zeros((txLUT['Naw'], nsinos), dtype=np.float32)
+    else:
+        assert output.shape == txLUT['Naw'], nsinos
+        assert output.dtype == np.dtype('float32')
     # fill the sino with gaps
-    mmr_auxe.rgaps(s, sino.astype(np.float32), txLUT, Cnt)
+    mmr_auxe.rgaps(output, cu.asarray(sino, dtype=np.float32), txLUT, Cnt)
 
     # return in the same data type as the input sino
-    return s.astype(sino.dtype)
+    return output.astype(sino.dtype)
 
 
 def mmrinit():
diff --git a/niftypet/nipet/src/aux_module.cu b/niftypet/nipet/src/aux_module.cu
index 3af8972d..e920d8f8 100644
--- a/niftypet/nipet/src/aux_module.cu
+++ b/niftypet/nipet/src/aux_module.cu
@@ -10,12 +10,13 @@ Copyrights: 2018
 #define PY_SSIZE_T_CLEAN
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION // NPY_API_VERSION
 
+#include "Python.h"
 #include "auxmath.h"
 #include "def.h"
 #include "norm.h"
+#include "numpy/arrayobject.h"
+#include "pycuvec.cuh"
 #include "scanner_0.h"
-#include <Python.h>
-#include <numpy/arrayobject.h>
 #include <stdlib.h>
 
 //=== START PYTHON INIT ===
@@ -386,26 +387,15 @@ static PyObject *mmr_pgaps(PyObject *self, PyObject *args) {
 
 //====================================================================================================
 static PyObject *mmr_rgaps(PyObject *self, PyObject *args) {
-
-  // output sino with gaps removed
-  PyObject *o_sng;
-
-  // transaxial LUT dictionary (e.g., 2D sino where dead bins are out).
-  PyObject *o_txLUT;
-
-  // Dictionary of scanner constants
-  PyObject *o_mmrcnst;
-
-  // input sino to be reformated with gaps removed
-  PyObject *o_sino;
-
-  // Structure of constants
-  Cnst Cnt;
-
-  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
-  /* Parse the input tuple */
-  if (!PyArg_ParseTuple(args, "OOOO", &o_sng, &o_sino, &o_txLUT, &o_mmrcnst)) return NULL;
-  //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+  PyCuVec<float> *o_sng = NULL;  // output sino with gaps removed
+  PyCuVec<float> *o_sino = NULL; // input sino to be reformated with gaps removed
+  PyObject *o_txLUT;   // transaxial LUT dictionary (e.g., 2D sino where dead bins are out).
+  PyObject *o_mmrcnst; // Dictionary of scanner constants
+  Cnst Cnt;            // Structure of constants
+  if (!PyArg_ParseTuple(args, "O&O&OO", &asPyCuVec_f, &o_sng, &asPyCuVec_f, &o_sino, &o_txLUT,
+                        &o_mmrcnst))
+    return NULL;
+  if (!o_sng || !o_sino) return NULL;
 
   /* Interpret the input objects as... PyLong_AsLong*/
   PyObject *pd_NSN11 = PyDict_GetItemString(o_mmrcnst, "NSN11");
@@ -427,45 +417,23 @@ static PyObject *mmr_rgaps(PyObject *self, PyObject *args) {
   PyObject *pd_aw2ali = PyDict_GetItemString(o_txLUT, "aw2ali");
 
   // input sino and the above 2D LUT
-  PyArrayObject *p_sino = NULL;
-  p_sino = (PyArrayObject *)PyArray_FROM_OTF(o_sino, NPY_FLOAT32, NPY_ARRAY_IN_ARRAY);
   PyArrayObject *p_aw2ali = NULL;
   p_aw2ali = (PyArrayObject *)PyArray_FROM_OTF(pd_aw2ali, NPY_INT32, NPY_ARRAY_IN_ARRAY);
+  if (p_aw2ali == NULL) { Py_XDECREF(p_aw2ali); }
+  int *aw2ali = (int *)PyArray_DATA(p_aw2ali);
 
   // number of sinogram from the shape of the sino (can be any number especially when using reduced
   // ring number)
-  int snno = (int)PyArray_DIM(p_sino, 0);
-
-  // output sino
-  PyArrayObject *p_sng = NULL;
-  p_sng = (PyArrayObject *)PyArray_FROM_OTF(o_sng, NPY_FLOAT32, NPY_ARRAY_INOUT_ARRAY2);
-
-  if (p_sino == NULL || p_aw2ali == NULL || p_sino == NULL) {
-    Py_XDECREF(p_aw2ali);
-    Py_XDECREF(p_sino);
-
-    PyArray_DiscardWritebackIfCopy(p_sng);
-    Py_XDECREF(p_sng);
-  }
-
-  int *aw2ali = (int *)PyArray_DATA(p_aw2ali);
-  float *sino = (float *)PyArray_DATA(p_sino);
-  float *sng = (float *)PyArray_DATA(p_sng);
-
-  // sets the device on which to calculate
-  HANDLE_ERROR(cudaSetDevice(Cnt.DEVID));
+  int snno = o_sino->shape[0];
 
   //<><><><><><><><><><><><><><><><><><><><><><>
+  HANDLE_ERROR(cudaSetDevice(Cnt.DEVID));
   // Run the conversion to GPU sinos
-  remove_gaps(sng, sino, snno, aw2ali, Cnt);
+  remove_gaps(o_sng->vec.data(), o_sino->vec.data(), snno, aw2ali, Cnt);
   //<><><><><><><><><><><><><><><><><><><><><><>
 
   // Clean up
   Py_DECREF(p_aw2ali);
-  Py_DECREF(p_sino);
-
-  PyArray_ResolveWritebackIfCopy(p_sng);
-  Py_DECREF(p_sng);
 
   Py_INCREF(Py_None);
   return Py_None;

From 97d9fbd0f7b0e7d9322bae8d4acc9e800dbb8f05 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sat, 27 Nov 2021 09:49:02 +0000
Subject: [PATCH 09/32] Revert "mmrprj: add `back_prj(div_sino=None)` CUDA
 convenience"

This reverts commit 1da5d5d39fc5457a35e3fffc199c3fa964c97266.
---
 niftypet/nipet/prj/mmrprj.py         | 11 ++-----
 niftypet/nipet/prj/src/prj_module.cu | 13 +++------
 niftypet/nipet/prj/src/prjb.cu       | 43 +++++++++-------------------
 niftypet/nipet/prj/src/prjb.h        |  2 +-
 4 files changed, 21 insertions(+), 48 deletions(-)

diff --git a/niftypet/nipet/prj/mmrprj.py b/niftypet/nipet/prj/mmrprj.py
index 2244d074..683155d4 100644
--- a/niftypet/nipet/prj/mmrprj.py
+++ b/niftypet/nipet/prj/mmrprj.py
@@ -144,8 +144,7 @@ def frwd_prj(im, scanner_params, isub=ISUB_DEFAULT, dev_out=False, attenuation=F
 # ------------------------------------------------------------------------
 
 
-def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=None, output=None,
-             sync=True):
+def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, output=None, sync=True):
     '''
     Calculate forward projection for the provided input image.
     Arguments:
@@ -154,7 +153,6 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=No
             transaxial and axial look up tables (LUT).
         isub -- array of transaxial indices of all sinograms (angles x bins) used for subsets;
             when the first element is negative, all transaxial bins are used (as in pure EM-ML).
-        div_sino -- if specificed, backprojects `sino[isub]/div_sino` instead of `sino`.
         output(CuVec, optional) -- output image.
         sync(bool) -- whether to `cudaDeviceSynchronize()` after.
     '''
@@ -180,10 +178,6 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=No
 
     # > check first the Siemens default sinogram;
     # > for this default shape only full sinograms are expected--no subsets.
-    orig_sino = sino
-    if div_sino is not None:
-        sino = sino[isub, :]
-        div_sino = cu.asarray(div_sino)
     if len(sino.shape) == 3:
         if sino.shape[0] != nsinos or sino.shape[1] != Cnt['NSANGLES'] or sino.shape[2] != Cnt[
                 'NSBINS']:
@@ -218,8 +212,7 @@ def back_prj(sino, scanner_params, isub=ISUB_DEFAULT, dev_out=False, div_sino=No
         assert bimg.dtype == np.dtype('float32')
 
     # > run back-projection
-    petprj.bprj(bimg, cu.asarray(sinog if div_sino is None else orig_sino), txLUT, axLUT, isub,
-                Cnt, div_sino=div_sino, sync=sync)
+    petprj.bprj(bimg, cu.asarray(sinog), txLUT, axLUT, isub, Cnt, sync=sync)
 
     if not dev_out:
         # > change from GPU optimised image dimensions to the standard Siemens shape
diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index ca05a0c8..c33dd464 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -424,18 +424,14 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   // output backprojected image
   PyCuVec<float> *o_bimg = NULL;
 
-  // sinogram divisor
-  PyCuVec<float> *o_div_sino = NULL;
-
   bool SYNC = true; // whether to ensure deviceToHost copy on return
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
   /* Parse the input tuple */
 
-  static const char *kwds[] = {"bimg", "sino",     "txLUT", "axLUT", "subs",
-                               "cnst", "div_sino", "sync",  NULL};
-  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&O&OOOO|O&b", (char **)kwds, &asPyCuVec_f,
+  static const char *kwds[] = {"bimg", "sino", "txLUT", "axLUT", "subs", "cnst", "sync", NULL};
+  if (!PyArg_ParseTupleAndKeywords(args, kwargs, "O&O&OOOO|b", (char **)kwds, &asPyCuVec_f,
                                    &o_bimg, &asPyCuVec_f, &o_sino, &o_txLUT, &o_axLUT, &o_subs,
-                                   &o_mmrcnst, &asPyCuVec_f, &o_div_sino, &SYNC))
+                                   &o_mmrcnst, &SYNC))
     return NULL;
   //^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
@@ -560,8 +556,7 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   HANDLE_ERROR(cudaMemcpy(d_subs, subs, Nprj * sizeof(int), cudaMemcpyHostToDevice));
 
   gpu_bprj(o_bimg->vec.data(), o_sino->vec.data(), li2rng, li2sn, li2nos, (short2 *)d_s2c,
-           (float4 *)d_crs, d_subs, d_tt, d_tv, Nprj, Cnt,
-           o_div_sino ? o_div_sino->vec.data() : nullptr, SYNC);
+           (float4 *)d_crs, d_subs, d_tt, d_tv, Nprj, Cnt, SYNC);
 
   HANDLE_ERROR(cudaFree(d_subs));
   HANDLE_ERROR(cudaFree(d_tv));
diff --git a/niftypet/nipet/prj/src/prjb.cu b/niftypet/nipet/prj/src/prjb.cu
index 2af74ca8..8ada3d4d 100644
--- a/niftypet/nipet/prj/src/prjb.cu
+++ b/niftypet/nipet/prj/src/prjb.cu
@@ -30,18 +30,12 @@ __global__ void imReduce(float *imr, float *im, int vz0, int nvz) {
 //===============================================================
 
 //**************** DIRECT ***********************************
-__global__ void bprj_drct(const float *sino, float *im, float *div_sino, const float *tt,
-                          const unsigned char *tv, const int *subs, const short snno) {
+__global__ void bprj_drct(const float *sino, float *im, const float *tt, const unsigned char *tv,
+                          const int *subs, const short snno) {
   int ixt = subs[blockIdx.x]; // transaxial indx
   int ixz = threadIdx.x;      // axial (z)
 
-  float bin;
-  if (div_sino) {
-    // using full sino with subset divisor (divisor may be smaller)
-    bin = sino[c_li2sn[ixz].x + ixt * snno] / div_sino[c_li2sn[ixz].x + blockIdx.x * snno];
-  } else {
-    bin = sino[c_li2sn[ixz].x + blockIdx.x * snno];
-  }
+  float bin = sino[c_li2sn[ixz].x + blockIdx.x * snno];
 
   float z = c_li2rng[ixz].x + .5 * SZ_RING;
   int w = (floorf(.5 * SZ_IMZ + SZ_VOXZi * z));
@@ -89,9 +83,8 @@ __global__ void bprj_drct(const float *sino, float *im, float *div_sino, const f
 }
 
 //************** OBLIQUE **************************************************
-__global__ void bprj_oblq(const float *sino, float *im, float *div_sino, const float *tt,
-                          const unsigned char *tv, const int *subs, const short snno,
-                          const int zoff, const short nil2r_c) {
+__global__ void bprj_oblq(const float *sino, float *im, const float *tt, const unsigned char *tv,
+                          const int *subs, const short snno, const int zoff, const short nil2r_c) {
 
   int ixz = threadIdx.x + zoff; // axial (z)
 
@@ -99,16 +92,8 @@ __global__ void bprj_oblq(const float *sino, float *im, float *div_sino, const f
 
     int ixt = subs[blockIdx.x]; // blockIdx.x is the transaxial bin index
                                 // bin values to be back projected
-
-    float bin, bin_;
-    if (div_sino) {
-      // using full sino with subset divisor (divisor may be smaller)
-      bin = sino[c_li2sn[ixz].x + snno * ixt] / div_sino[c_li2sn[ixz].x + snno * blockIdx.x];
-      bin_ = sino[c_li2sn[ixz].y + snno * ixt] / div_sino[c_li2sn[ixz].y + snno * blockIdx.x];
-    } else {
-      bin = sino[c_li2sn[ixz].x + snno * blockIdx.x];
-      bin_ = sino[c_li2sn[ixz].y + snno * blockIdx.x];
-    }
+    float bin = sino[c_li2sn[ixz].x + snno * blockIdx.x];
+    float bin_ = sino[c_li2sn[ixz].y + snno * blockIdx.x];
 
     //-------------------------------------------------
     /*** accumulation ***/
@@ -204,7 +189,7 @@ __global__ void bprj_oblq(const float *sino, float *im, float *div_sino, const f
 //--------------------------------------------------------------------------------------------------
 void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short2 *d_s2c,
               float4 *d_crs, int *d_subs, float *d_tt, unsigned char *d_tv, int Nprj, Cnst Cnt,
-              float *_d_div_sino, bool _sync) {
+              bool _sync) {
   int dev_id;
   cudaGetDevice(&dev_id);
   if (Cnt.LOG <= LOGDEBUG) printf("i> using CUDA device #%d\n", dev_id);
@@ -272,7 +257,7 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
   //-----------------------------------------------------------------------
 
   //============================================================================
-  bprj_drct<<<Nprj, nrng_c>>>(d_sino, d_imf, _d_div_sino, d_tt, d_tv, d_subs, snno);
+  bprj_drct<<<Nprj, nrng_c>>>(d_sino, d_imf, d_tt, d_tv, d_subs, snno);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
@@ -282,11 +267,11 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
   int Nz = ((Noblq + 127) / 128) * 128;
 
   //============================================================================
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, _d_div_sino, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
   HANDLE_ERROR(cudaGetLastError());
 
   zoff += Nz / 2;
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, _d_div_sino, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_imf, d_tt, d_tv, d_subs, snno, zoff, nil2r_c);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
@@ -356,19 +341,19 @@ void rec_bprj(float *d_bimg, float *d_sino, int *d_sub, int Nprj, float *d_tt, u
   if (Cnt.LOG <= LOGDEBUG) printf("i> subset back projection (Nprj=%d)... ", Nprj);
 
   //============================================================================
-  bprj_drct<<<Nprj, NRINGS>>>(d_sino, d_bimg, nullptr, d_tt, d_tv, d_sub, snno);
+  bprj_drct<<<Nprj, NRINGS>>>(d_sino, d_bimg, d_tt, d_tv, d_sub, snno);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
   int zoff = NRINGS;
   //============================================================================
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, nullptr, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
   zoff += Nz / 2;
   //============================================================================
-  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, nullptr, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
+  bprj_oblq<<<Nprj, Nz / 2>>>(d_sino, d_bimg, d_tt, d_tv, d_sub, snno, zoff, NLI2R);
   HANDLE_ERROR(cudaGetLastError());
   //============================================================================
 
diff --git a/niftypet/nipet/prj/src/prjb.h b/niftypet/nipet/prj/src/prjb.h
index 447392bb..0237e234 100644
--- a/niftypet/nipet/prj/src/prjb.h
+++ b/niftypet/nipet/prj/src/prjb.h
@@ -9,7 +9,7 @@
 // used from Python
 void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2nos, short2 *d_s2c,
               float4 *d_crs, int *d_subs, float *d_tt, unsigned char *d_tv, int Nprj, Cnst Cnt,
-              float *_d_div_sino = nullptr, bool _sync = true);
+              bool _sync = true);
 
 // to be used within CUDA C reconstruction
 void rec_bprj(float *d_bimg, float *d_sino, int *sub, int Nprj,

From 23338ecfa946e53568d307de39264a9d59f504c1 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sat, 27 Nov 2021 18:19:21 +0000
Subject: [PATCH 10/32] remove unnecessary NULL checks

---
 niftypet/nipet/prj/src/prj_module.cu | 5 ++---
 niftypet/nipet/src/aux_module.cu     | 1 -
 2 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index c33dd464..28a6eb6b 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -315,8 +315,7 @@ static PyObject *frwd_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
   /* If that didn't work, throw an exception. */
   if (p_li2rno == NULL || p_li2sn == NULL || p_li2sn1 == NULL || p_li2nos == NULL ||
-      p_aw2ali == NULL || p_s2c == NULL || !o_im || p_crs == NULL || p_subs == NULL || !o_prjout ||
-      p_li2rng == NULL) {
+      p_aw2ali == NULL || p_s2c == NULL || p_crs == NULL || p_subs == NULL || p_li2rng == NULL) {
     // axLUTs
     Py_XDECREF(p_li2rno);
     Py_XDECREF(p_li2sn);
@@ -480,7 +479,7 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
 
   /* If that didn't work, throw an exception. */
   if (p_li2rno == NULL || p_li2sn == NULL || p_li2sn1 == NULL || p_li2nos == NULL ||
-      p_s2c == NULL || !o_sino || p_crs == NULL || p_subs == NULL || p_li2rng == NULL || !o_bimg) {
+      p_s2c == NULL || p_crs == NULL || p_subs == NULL || p_li2rng == NULL) {
     // axLUTs
     Py_XDECREF(p_li2rno);
     Py_XDECREF(p_li2sn);
diff --git a/niftypet/nipet/src/aux_module.cu b/niftypet/nipet/src/aux_module.cu
index e920d8f8..d6be873b 100644
--- a/niftypet/nipet/src/aux_module.cu
+++ b/niftypet/nipet/src/aux_module.cu
@@ -395,7 +395,6 @@ static PyObject *mmr_rgaps(PyObject *self, PyObject *args) {
   if (!PyArg_ParseTuple(args, "O&O&OO", &asPyCuVec_f, &o_sng, &asPyCuVec_f, &o_sino, &o_txLUT,
                         &o_mmrcnst))
     return NULL;
-  if (!o_sng || !o_sino) return NULL;
 
   /* Interpret the input objects as... PyLong_AsLong*/
   PyObject *pd_NSN11 = PyDict_GetItemString(o_mmrcnst, "NSN11");

From c281ba780f7c8bb7d4dcca52f2acba2a18e6a28a Mon Sep 17 00:00:00 2001
From: Pawel <p.markiewicz@ucl.ac.uk>
Date: Thu, 9 Dec 2021 01:14:08 +0000
Subject: [PATCH 11/32] enhanced trimming withing rocnsturction when using
 pipe.py

---
 niftypet/nipet/img/pipe.py | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/niftypet/nipet/img/pipe.py b/niftypet/nipet/img/pipe.py
index df94e994..8c34ab41 100644
--- a/niftypet/nipet/img/pipe.py
+++ b/niftypet/nipet/img/pipe.py
@@ -47,6 +47,8 @@ def mmrchain(
     trim_interp=0,          # interpolation for upsampling used in PVC
     trim_memlim=True,       # reduced use of memory for machines
                             # with limited memory (slow though)
+    trim_store=True,
+
     pvcroi=None,            # ROI used for PVC.  If undefined no PVC
                             # is performed.
     pvcreg_tool='niftyreg', # the registration tool used in PVC
@@ -406,8 +408,15 @@ def mmrchain(
         elif 'lm_ima' in datain:
             fnm = os.path.basename(datain['lm_ima'])[:20]
         # trim PET and upsample
-        petu = nimpa.imtrimup(dynim, affine=image_affine(datain, Cnt), scale=trim_scale,
-                              int_order=trim_interp, outpath=petimg, fname=fnm, fcomment=fcomment,
+        petu = nimpa.imtrimup(dynim,
+                              affine=image_affine(datain, Cnt),
+                              scale=trim_scale,
+                              int_order=trim_interp,
+                              outpath=petimg,
+                              fname=fnm,
+                              fcomment=fcomment,
+                              fcomment_pfx=fout+'_',
+                              store_img=trim_store, 
                               store_img_intrmd=store_img_intrmd, memlim=trim_memlim,
                               verbose=log.getEffectiveLevel())
 

From c5df379487daadeefb4203891cde2c64df6576b4 Mon Sep 17 00:00:00 2001
From: Pawel <p.markiewicz@ucl.ac.uk>
Date: Thu, 9 Dec 2021 01:28:13 +0000
Subject: [PATCH 12/32] fixing storing trimming naming when reconstructing

---
 niftypet/nipet/img/pipe.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/niftypet/nipet/img/pipe.py b/niftypet/nipet/img/pipe.py
index 8c34ab41..5ca1b8dd 100644
--- a/niftypet/nipet/img/pipe.py
+++ b/niftypet/nipet/img/pipe.py
@@ -548,11 +548,11 @@ def mmrchain(
             if fout is None:
                 fpetu = os.path.join(
                     pettrim,
-                    os.path.basename(fpet) + f'_trimmed-upsampled-scale-{trim_scale}')
+                    os.path.basename(fpet) + f'_recon-trimmed-upsampled-scale-{trim_scale}')
             else:
                 fpetu = os.path.join(
                     pettrim,
-                    os.path.basename(fout) + f'_trimmed-upsampled-scale-{trim_scale}')
+                    os.path.basename(fout) + f'_recon-trimmed-upsampled-scale-{trim_scale}')
             # in case of PVC
             if pvcroi:
                 # itertive Yang (iY) added to NIfTI descritoption

From fc0923465ff3ae8ef08ded77b4e0983e4ccc8fa3 Mon Sep 17 00:00:00 2001
From: Pawel <p.markiewicz@ucl.ac.uk>
Date: Sat, 25 Dec 2021 00:12:17 +0000
Subject: [PATCH 13/32] type in comments

---
 niftypet/nipet/img/mmrimg.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index cff4e6c0..6ccefa9d 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -1101,7 +1101,7 @@ def get_hmupos(datain, parts, Cnt, outpath=''):
     nimpa.array2nii(np.zeros((Cnt['SO_IMZ'], Cnt['SO_IMY'], Cnt['SO_IMX']), dtype=np.float32), B,
                     fref)
 
-    # pefine a dictionary of all positions/offsets of hardware mu-maps
+    # define a dictionary of all positions/offsets of hardware mu-maps
     hmupos = [None] * 5
     hmupos[0] = {
         'TabPosOrg': tpozyx, # prom DICOM of LM file

From 35d6bdfb20ffc68c0c90cc9b06b08857dab0116d Mon Sep 17 00:00:00 2001
From: Pawel <p.markiewicz@ucl.ac.uk>
Date: Fri, 31 Dec 2021 14:34:32 +0000
Subject: [PATCH 14/32] changed resampling method of the mu-maps

---
 niftypet/nipet/img/mmrimg.py | 53 ++++++++++++++++--------------------
 niftypet/nipet/prj/mmrrec.py |  3 +-
 2 files changed, 25 insertions(+), 31 deletions(-)

diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index 6ccefa9d..4b755cc7 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -279,15 +279,13 @@ def mudcm2nii(datain, Cnt):
     nimpa.array2nii(im, B, os.path.join(os.path.dirname(datain['mumapDCM']), 'muref.nii.gz'))
     # -------------------------------------------------------------------------------------
     fmu = os.path.join(os.path.dirname(datain['mumapDCM']), 'mu_r.nii.gz')
-    if os.path.isfile(Cnt['RESPATH']):
-        run([
-            Cnt['RESPATH'], '-ref',
-            os.path.join(os.path.dirname(datain['mumapDCM']), 'muref.nii.gz'), '-flo',
-            os.path.join(os.path.dirname(datain['mumapDCM']), 'mu.nii.gz'), '-res', fmu, '-pad',
-            '0'])
-    else:
-        log.error('path to resampling executable is incorrect!')
-        raise IOError('Error launching NiftyReg for image resampling.')
+
+    nimpa.resample_dipy(
+        os.path.join(os.path.dirname(datain['mumapDCM']), 'muref.nii.gz'),
+        os.path.join(os.path.dirname(datain['mumapDCM']), 'mu.nii.gz'),
+        fimout=fmu,
+        intrp=1,
+        dtype_nifti=np.float32)
 
     return fmu
 
@@ -346,23 +344,21 @@ def obj_mumap(
     for d in resdcm:
         os.remove(d)
 
-    # convert the DICOM mu-map images to nii
+    # > convert the DICOM mu-map images to NIfTI
     run([Cnt['DCM2NIIX'], '-f', fnii + tstmp, '-o', fmudir, datain['mumapDCM']])
     # piles for the T1w, pick one:
     fmunii = glob.glob(os.path.join(fmudir, '*' + fnii + tstmp + '*.nii*'))[0]
     # fmunii = glob.glob( os.path.join(datain['mumapDCM'], '*converted*.nii*') )
     # fmunii = fmunii[0]
 
-    # the converted nii image resample to the reference size
+    # > resampled the NIfTI converted image to the reference shape/size
     fmu = os.path.join(fmudir, comment + 'mumap_tmp.nii.gz')
-    if os.path.isfile(Cnt['RESPATH']):
-        cmd = [Cnt['RESPATH'], '-ref', fmuref, '-flo', fmunii, '-res', fmu, '-pad', '0']
-        if log.getEffectiveLevel() > logging.INFO:
-            cmd.append('-voff')
-        run(cmd)
-    else:
-        log.error('path to resampling executable is incorrect!')
-        raise IOError('Path to executable is incorrect!')
+    nimpa.resample_dipy(
+        fmuref,
+        fmunii,
+        fimout=fmu,
+        intrp=1,
+        dtype_nifti=np.float32)
 
     nim = nib.load(fmu)
     # get the affine transform
@@ -1034,9 +1030,7 @@ def rd_hmu(fh):
 
 
 def get_hmupos(datain, parts, Cnt, outpath=''):
-    # check if registration executable exists
-    if not os.path.isfile(Cnt['RESPATH']):
-        raise IOError('No registration executable found!')
+
     # ----- get positions from the DICOM list-mode file -----
     ihdr, csainfo = mmraux.hdr_lm(datain, Cnt)
     # pable position origin
@@ -1142,15 +1136,16 @@ def get_hmupos(datain, parts, Cnt, outpath=''):
         A[2, 3] = -10 * ((s[0] - 1) * vs[0] - vpos[0])
         nimpa.array2nii(im[::-1, ::-1, :], A, hmupos[i]['niipath'])
 
-        # resample using nify.reg
+        # > resample using DIPY in nimpa
         fout = os.path.join(os.path.dirname(hmupos[0]['niipath']),
                             'r' + os.path.basename(hmupos[i]['niipath']).split('.')[0] + '.nii.gz')
-        cmd = [
-            Cnt['RESPATH'], '-ref', hmupos[0]['niipath'], '-flo', hmupos[i]['niipath'], '-res',
-            fout, '-pad', '0']
-        if log.getEffectiveLevel() > logging.INFO:
-            cmd.append('-voff')
-        run(cmd)
+
+        nimpa.resample_dipy(
+            hmupos[0]['niipath'],
+            hmupos[i]['niipath'],
+            fimout=fout,
+            intrp=1,
+            dtype_nifti=np.float32)
 
     return hmupos
 
diff --git a/niftypet/nipet/prj/mmrrec.py b/niftypet/nipet/prj/mmrrec.py
index a2b2411d..841469cd 100644
--- a/niftypet/nipet/prj/mmrrec.py
+++ b/niftypet/nipet/prj/mmrrec.py
@@ -232,8 +232,7 @@ def osemone(datain, mumaps, hst, scanner_params, recmod=3, itr=4, fwhm=0., psf=N
             log.info('using provided attenuation factor sinogram')
         else:
             asng = cu.zeros(psng.shape, dtype=np.float32)
-            petprj.fprj(asng, cu.asarray(mus), txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt,
-                        1)
+            petprj.fprj(asng, cu.asarray(mus), txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt, 1)
     # > combine attenuation and normalisation
     ansng = asng * nsng
     # ========================================================================

From 733c039bd31ba621776e81b1e68fee647fe0e82a Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sun, 28 Nov 2021 08:59:08 +0000
Subject: [PATCH 15/32] bubble up prjf mallocs

---
 niftypet/nipet/prj/src/prj_module.cu | 35 +++++++++++++++++++---
 niftypet/nipet/prj/src/prjb.h        |  9 ++----
 niftypet/nipet/prj/src/prjf.cu       | 43 ++--------------------------
 niftypet/nipet/prj/src/prjf.h        |  6 ++--
 4 files changed, 39 insertions(+), 54 deletions(-)

diff --git a/niftypet/nipet/prj/src/prj_module.cu b/niftypet/nipet/prj/src/prj_module.cu
index 28a6eb6b..b54e7362 100644
--- a/niftypet/nipet/prj/src/prj_module.cu
+++ b/niftypet/nipet/prj/src/prj_module.cu
@@ -376,8 +376,35 @@ static PyObject *frwd_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   HANDLE_ERROR(cudaSetDevice(Cnt.DEVID));
 
   //<><><><><><><<><><><><><><><><><><><><><><><><><<><><><><><><><><><><><><><><><><><><<><><><><><><><><><><>
-  gpu_fprj(o_prjout->vec.data(), o_im->vec.data(), li2rng, li2sn, li2nos, s2c, aw2ali, crs, subs,
-           Nprj, Naw, N0crs, Cnt, att, SYNC);
+  //--- TRANSAXIAL COMPONENT
+  float4 *d_crs;
+  HANDLE_ERROR(cudaMalloc(&d_crs, N0crs * sizeof(float4)));
+  HANDLE_ERROR(cudaMemcpy(d_crs, crs, N0crs * sizeof(float4), cudaMemcpyHostToDevice));
+
+  short2 *d_s2c;
+  HANDLE_ERROR(cudaMalloc(&d_s2c, AW * sizeof(short2)));
+  HANDLE_ERROR(cudaMemcpy(d_s2c, s2c, AW * sizeof(short2), cudaMemcpyHostToDevice));
+
+  float *d_tt;
+  HANDLE_ERROR(cudaMalloc(&d_tt, N_TT * AW * sizeof(float)));
+
+  unsigned char *d_tv;
+  HANDLE_ERROR(cudaMalloc(&d_tv, N_TV * AW * sizeof(unsigned char)));
+  HANDLE_ERROR(cudaMemset(d_tv, 0, N_TV * AW * sizeof(unsigned char)));
+
+  // array of subset projection bins
+  int *d_subs;
+  HANDLE_ERROR(cudaMalloc(&d_subs, Nprj * sizeof(int)));
+  HANDLE_ERROR(cudaMemcpy(d_subs, subs, Nprj * sizeof(int), cudaMemcpyHostToDevice));
+
+  gpu_fprj(o_prjout->vec.data(), o_im->vec.data(), li2rng, li2sn, li2nos, d_s2c, aw2ali, d_crs,
+           d_subs, d_tt, d_tv, Nprj, Naw, Cnt, att, SYNC);
+
+  HANDLE_ERROR(cudaFree(d_subs));
+  HANDLE_ERROR(cudaFree(d_tv));
+  HANDLE_ERROR(cudaFree(d_tt));
+  HANDLE_ERROR(cudaFree(d_s2c));
+  HANDLE_ERROR(cudaFree(d_crs));
   //<><><><><><><><<><><><><><><><><><><><><><><><><<><><><><><><><><><><><><><><><><><><<><><><><><><><><><><>
 
   // Clean up
@@ -554,8 +581,8 @@ static PyObject *back_prj(PyObject *self, PyObject *args, PyObject *kwargs) {
   HANDLE_ERROR(cudaMalloc(&d_subs, Nprj * sizeof(int)));
   HANDLE_ERROR(cudaMemcpy(d_subs, subs, Nprj * sizeof(int), cudaMemcpyHostToDevice));
 
-  gpu_bprj(o_bimg->vec.data(), o_sino->vec.data(), li2rng, li2sn, li2nos, (short2 *)d_s2c,
-           (float4 *)d_crs, d_subs, d_tt, d_tv, Nprj, Cnt, SYNC);
+  gpu_bprj(o_bimg->vec.data(), o_sino->vec.data(), li2rng, li2sn, li2nos, d_s2c, d_crs, d_subs,
+           d_tt, d_tv, Nprj, Cnt, SYNC);
 
   HANDLE_ERROR(cudaFree(d_subs));
   HANDLE_ERROR(cudaFree(d_tv));
diff --git a/niftypet/nipet/prj/src/prjb.h b/niftypet/nipet/prj/src/prjb.h
index 0237e234..0c018965 100644
--- a/niftypet/nipet/prj/src/prjb.h
+++ b/niftypet/nipet/prj/src/prjb.h
@@ -12,12 +12,7 @@ void gpu_bprj(float *d_im, float *d_sino, float *li2rng, short *li2sn, char *li2
               bool _sync = true);
 
 // to be used within CUDA C reconstruction
-void rec_bprj(float *d_bimg, float *d_sino, int *sub, int Nprj,
-
-              float *d_tt, unsigned char *d_tv,
-
-              float *li2rng, short *li2sn, char *li2nos,
-
-              Cnst Cnt);
+void rec_bprj(float *d_bimg, float *d_sino, int *sub, int Nprj, float *d_tt, unsigned char *d_tv,
+              float *li2rng, short *li2sn, char *li2nos, Cnst Cnt);
 
 #endif
diff --git a/niftypet/nipet/prj/src/prjf.cu b/niftypet/nipet/prj/src/prjf.cu
index 9c937c79..5fe70325 100644
--- a/niftypet/nipet/prj/src/prjf.cu
+++ b/niftypet/nipet/prj/src/prjf.cu
@@ -206,35 +206,13 @@ __global__ void fprj_oblq(float *sino, const float *im, const float *tt, const u
 }
 
 //--------------------------------------------------------------------------------------------------
-void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt, char att,
-              bool _sync) {
+void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2nos, short2 *d_s2c,
+              int *aw2ali, float4 *d_crs, int *d_subs, float *d_tt, unsigned char *d_tv, int Nprj,
+              int Naw, Cnst Cnt, char att, bool _sync) {
   int dev_id;
   cudaGetDevice(&dev_id);
   if (Cnt.LOG <= LOGDEBUG) printf("i> using CUDA device #%d\n", dev_id);
 
-  //--- TRANSAXIAL COMPONENT
-  float4 *d_crs;
-  HANDLE_ERROR(cudaMalloc(&d_crs, N0crs * sizeof(float4)));
-  HANDLE_ERROR(cudaMemcpy(d_crs, crs, N0crs * sizeof(float4), cudaMemcpyHostToDevice));
-
-  short2 *d_s2c;
-  HANDLE_ERROR(cudaMalloc(&d_s2c, AW * sizeof(short2)));
-  HANDLE_ERROR(cudaMemcpy(d_s2c, s2c, AW * sizeof(short2), cudaMemcpyHostToDevice));
-
-  float *d_tt;
-  HANDLE_ERROR(cudaMalloc(&d_tt, N_TT * AW * sizeof(float)));
-
-  unsigned char *d_tv;
-  HANDLE_ERROR(cudaMalloc(&d_tv, N_TV * AW * sizeof(unsigned char)));
-  HANDLE_ERROR(cudaMemset(d_tv, 0, N_TV * AW * sizeof(unsigned char)));
-
-  // array of subset projection bins
-  int *d_subs;
-  HANDLE_ERROR(cudaMalloc(&d_subs, Nprj * sizeof(int)));
-  HANDLE_ERROR(cudaMemcpy(d_subs, subs, Nprj * sizeof(int), cudaMemcpyHostToDevice));
-  //---
-
   //-----------------------------------------------------------------
   // RINGS: either all or a subset of rings can be used (span-1 feature only)
   //-----------------------------------------------------------------
@@ -289,16 +267,6 @@ void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2no
     HANDLE_ERROR(cudaGetLastError());
   }
 
-  // float *d_li2rng;  HANDLE_ERROR( cudaMalloc(&d_li2rng, N0li*N1li*sizeof(float)) );
-  // HANDLE_ERROR( cudaMemcpy( d_li2rng, li2rng, N0li*N1li*sizeof(float), cudaMemcpyHostToDevice)
-  // );
-
-  // int *d_li2sn;  HANDLE_ERROR(cudaMalloc(&d_li2sn, N0li*N1li*sizeof(int)) );
-  // HANDLE_ERROR( cudaMemcpy( d_li2sn, li2sn, N0li*N1li*sizeof(int), cudaMemcpyHostToDevice) );
-
-  // int *d_li2nos;  HANDLE_ERROR( cudaMalloc(&d_li2nos, N1li*sizeof(int)) );
-  // HANDLE_ERROR( cudaMemcpy( d_li2nos, li2nos, N1li*sizeof(int), cudaMemcpyHostToDevice) );
-
   cudaMemcpyToSymbol(c_li2rng, li2rng, nil2r_c * sizeof(float2));
   cudaMemcpyToSymbol(c_li2sn, li2sn, nil2r_c * sizeof(short2));
   cudaMemcpyToSymbol(c_li2nos, li2nos, nil2r_c * sizeof(char));
@@ -349,11 +317,6 @@ void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2no
   }
 
   if (nvz < SZ_IMZ) HANDLE_ERROR(cudaFree(d_im));
-  HANDLE_ERROR(cudaFree(d_tt));
-  HANDLE_ERROR(cudaFree(d_tv));
-  HANDLE_ERROR(cudaFree(d_subs));
-  HANDLE_ERROR(cudaFree(d_crs));
-  HANDLE_ERROR(cudaFree(d_s2c));
 }
 
 //=======================================================================
diff --git a/niftypet/nipet/prj/src/prjf.h b/niftypet/nipet/prj/src/prjf.h
index 38bd525d..08a66294 100644
--- a/niftypet/nipet/prj/src/prjf.h
+++ b/niftypet/nipet/prj/src/prjf.h
@@ -6,9 +6,9 @@
 #ifndef PRJF_H
 #define PRJF_H
 
-void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2nos, short *s2c,
-              int *aw2ali, float *crs, int *subs, int Nprj, int Naw, int N0crs, Cnst Cnt, char att,
-              bool _sync = true);
+void gpu_fprj(float *d_sn, float *d_im, float *li2rng, short *li2sn, char *li2nos, short2 *d_s2c,
+              int *aw2ali, float4 *d_crs, int *d_subs, float *d_tt, unsigned char *d_tv, int Nprj,
+              int Naw, Cnst Cnt, char att, bool _sync = true);
 
 void rec_fprj(float *d_sino, float *d_img, int *d_sub, int Nprj,
 

From 783649656135b07fbf5c79e9d0a3a8f1818c3ef0 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Thu, 23 Dec 2021 12:52:56 +0000
Subject: [PATCH 16/32] tests: fix deprecation warning

---
 tests/test_amyloid_pvc.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/test_amyloid_pvc.py b/tests/test_amyloid_pvc.py
index a7bd979f..671a7986 100644
--- a/tests/test_amyloid_pvc.py
+++ b/tests/test_amyloid_pvc.py
@@ -66,7 +66,7 @@ def muhdct(mMRpars, datain, tmp_path_factory, worker_id):
 
     opth = str(tmp_path.parent / "muhdct")
     flock = FileLock(opth + ".lock")
-    with flock.acquire(poll_intervall=0.5): # xdist, force auto-reuse via flock
+    with flock.acquire(poll_interval=0.5): # xdist, force auto-reuse via flock
         return nipet.hdw_mumap(datain, [1, 2, 4], mMRpars, outpath=opth, use_stored=True)
 
 

From fac7983b5656b7ed36e85af9964b018743b85f78 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Thu, 23 Dec 2021 12:56:07 +0000
Subject: [PATCH 17/32] build: update pre-commit

---
 .pre-commit-config.yaml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index b5956c1f..77e4b347 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -41,7 +41,7 @@ repos:
   - id: yapf
     args: [-i]
 - repo: https://github.com/PyCQA/isort
-  rev: 5.9.3
+  rev: 5.10.1
   hooks:
   - id: isort
 - repo: https://github.com/doublify/pre-commit-clang-format

From f46f9f1469303911405d43abde93a8a528f5bdc7 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sun, 2 Jan 2022 11:32:21 +0000
Subject: [PATCH 18/32] examples: update demo notebook

---
 examples/demo.ipynb | 322 +++++++++++++++++++++++++++++++-------------
 setup.cfg           |   2 +-
 2 files changed, 227 insertions(+), 97 deletions(-)

diff --git a/examples/demo.ipynb b/examples/demo.ipynb
index afde8439..a03def75 100644
--- a/examples/demo.ipynb
+++ b/examples/demo.ipynb
@@ -12,10 +12,10 @@
     "Mathematically:\n",
     "\n",
     "$$\n",
-    "{\\bf\\theta}^{(k+1)} = {{\\bf\\theta}^{(k)} \\over \\sum_n{{\\bf H}^T{\\bf X}_n^T{\\bf A}^T{\\bf N}^T{\\bf 1}}}\n",
+    "{\\bf y}^{(k+1)} = {{\\bf y}^{(k)} \\over \\sum_n{{\\bf H}^T{\\bf X}_n^T{\\bf A}^T{\\bf N}^T{\\bf 1}}}\n",
     "    \\circ\n",
     "    \\sum_n{ {\\bf H}^T{\\bf X}_n^T{\\bf A}^T{\\bf N}^T\n",
-    "        { {\\bf m} \\over {\\bf NA}{\\bf X}_n{\\bf H}{\\bf\\theta}^{(k)} + {\\bf r} + {\\bf s} }\n",
+    "        { {\\bf m} \\over {\\bf NA}{\\bf X}_n{\\bf H}{\\bf y}^{(k)} + {\\bf r} + {\\bf s} }\n",
     "    },\n",
     "$$\n",
     "\n",
@@ -26,8 +26,8 @@
     "\n",
     "----\n",
     "\n",
-    "- Author: Casper O. da Costa-Luis [casper.dcl@{physics.org|ieee.org|kcl.ac.uk}](mailto:casper.dcl@physics.org)\n",
-    "- Date: 2019-20\n",
+    "- Author: Casper O. da Costa-Luis [casper.dcl@{physics.org|ieee.org|ucl.ac.uk|kcl.ac.uk}](mailto:casper.dcl@physics.org)\n",
+    "- Date: 2019-21\n",
     "\n",
     "----\n",
     "\n",
@@ -41,29 +41,30 @@
    "outputs": [],
    "source": [
     "# imports\n",
-    "from __future__ import print_function, division\n",
     "%matplotlib notebook\n",
     "\n",
-    "from niftypet import nipet\n",
-    "from niftypet import nimpa\n",
+    "import logging\n",
+    "from functools import partial\n",
+    "from os import path\n",
+    "from pathlib import Path\n",
     "\n",
+    "import cuvec as cu\n",
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
     "from scipy.ndimage.filters import gaussian_filter\n",
     "from tqdm.auto import trange\n",
-    "from brainweb import volshow\n",
     "\n",
-    "from collections import OrderedDict\n",
-    "from os import path\n",
-    "import functools\n",
-    "import logging\n",
-    "logging.basicConfig(level=logging.INFO, format=nipet.LOG_FORMAT)\n",
+    "from niftypet import nipet, nimpa\n",
     "\n",
+    "logging.basicConfig(level=logging.INFO, format=nipet.LOG_FORMAT)\n",
     "print(nipet.gpuinfo())\n",
     "# get all the scanner parameters\n",
     "mMRpars = nipet.get_mmrparams()\n",
+    "mMRpars['Cnt']['LOG'] = logging.INFO\n",
     "# conversion for Gaussian sigma/[voxel] to FWHM/[mm]\n",
-    "SIGMA2FWHMmm = (8 * np.log(2))**0.5 * np.array([mMRpars['Cnt']['SO_VX' + i] for i in 'ZYX']) * 10"
+    "SIGMA2FWHMmm = (8 * np.log(2))**0.5 * np.array([mMRpars['Cnt']['SO_VX' + i] for i in 'ZYX']) * 10\n",
+    "# image-space PSF function\n",
+    "imsmooth = partial(nimpa.imsmooth, Cnt=mMRpars['Cnt'])"
    ]
   },
   {
@@ -79,17 +80,12 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "folderin = \"amyloidPET_FBP_TP0\"\n",
-    "\n",
+    "folderin = Path(\"amyloidPET_FBP_TP0\")\n",
     "# automatically categorise the input data\n",
-    "#logging.getLogger().setLevel(logging.INFO)\n",
     "datain = nipet.classify_input(folderin, mMRpars)\n",
-    "\n",
     "# output path\n",
-    "opth = path.join(datain['corepath'], 'niftyout')\n",
-    "# switch on verbose mode\n",
-    "#logging.getLogger().setLevel(logging.DEBUG)\n",
-    "\n",
+    "opth = Path(datain['corepath']) / \"niftyout\"\n",
+    "# show categorised data\n",
     "datain"
    ]
   },
@@ -100,7 +96,23 @@
    "outputs": [],
    "source": [
     "# hardware mu-map (bed, head/neck coils)\n",
-    "mu_h = nipet.hdw_mumap(datain, [1,2,4], mMRpars, outpath=opth, use_stored=True)"
+    "mu_h = nipet.hdw_mumap(datain, [1, 2, 4], mMRpars, outpath=opth, use_stored=True)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# create histogram\n",
+    "mMRpars['Cnt']['BTP'] = 0\n",
+    "mSino = nipet.mmrhist(datain, mMRpars, outpath=opth, store=True, use_stored=True)\n",
+    "if False:                                                   # enable parametric bootstrap\n",
+    "    mMRpars['Cnt']['BTP'] = 2\n",
+    "    totCnt = 3e6                                            # total counts to simulate\n",
+    "    mMRpars['Cnt']['BTPRT'] = totCnt / mSino['psino'].sum() # count level ratio relative to original\n",
+    "    mSino = nipet.mmrhist(datain, mMRpars, outpath=opth / \"BTP\" / f\"{totCnt:.3g}\", store=True)"
    ]
   },
   {
@@ -110,22 +122,23 @@
    "outputs": [],
    "source": [
     "# MR-based human mu-map\n",
-    "\n",
-    "# UTE-based object mu-map aligned (need UTE sequence or T1 for pseudo-CT)\n",
-    "#mu_o = nipet.align_mumap(\n",
-    "#    datain,\n",
-    "#    scanner_params=mMRpars,\n",
-    "#    outpath=opth,\n",
-    "#    t0=0, t1=0, # when both times are 0, will use full data\n",
-    "#    itr=2,      # number of iterations used for recon to which registering MR/UTE\n",
-    "#    petopt='ac',# what PET image to use (ac-just attenuation corrected)\n",
-    "#    musrc='ute',# source of mu-map (ute/pct)\n",
-    "#    ute_name='UTE2', # which UTE to use (UTE1/2 shorter/faster)\n",
-    "#    verbose=True,\n",
-    "#)\n",
-    "\n",
-    "#> the same as above without any faff though (no alignment)\n",
-    "mu_o = nipet.obj_mumap(datain, mMRpars, outpath=opth, store=True)"
+    "if True:                        # UTE-based object mu-map aligned (need UTE sequence or T1 for pseudo-CT)\n",
+    "    mu_o = nipet.align_mumap(\n",
+    "        datain,\n",
+    "        scanner_params=mMRpars,\n",
+    "        outpath=opth,\n",
+    "        store=True,\n",
+    "        use_stored=True,\n",
+    "        hst=mSino,\n",
+    "        t0=0,\n",
+    "        t1=0,                   # when both times are 0, will use full data\n",
+    "        itr=2,                  # number of iterations used for recon to which registering MR/UTE\n",
+    "        petopt='ac',            # what PET image to use (ac-just attenuation corrected)\n",
+    "        musrc='pct',            # source of mu-map (ute/pct)\n",
+    "        ute_name='UTE2',        # which UTE to use (UTE1/2 shorter/faster)\n",
+    "    )\n",
+    "else:                           # the same as above without any faff though (no alignment)\n",
+    "    mu_o = nipet.obj_mumap(datain, mMRpars, outpath=opth, store=True)"
    ]
   },
   {
@@ -134,14 +147,11 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# create histogram\n",
-    "mMRpars['Cnt']['BTP'] = 0\n",
-    "m = nipet.mmrhist(datain, mMRpars, outpath=opth, store=True, use_stored=True)\n",
-    "if False:\n",
-    "    mMRpars['Cnt']['BTP'] = 2  # enable parametric bootstrap\n",
-    "    totCnt = 3e6\n",
-    "    mMRpars['Cnt']['BTPRT'] = totCnt / m['psino'].sum()  # ratio count level relative to the original\n",
-    "    m = nipet.mmrhist(datain, mMRpars, outpath=path.join(opth, 'BTP', '%.3g' % totCnt), store=True)"
+    "try:\n",
+    "    mu = mu_h['im'] + mu_o['im'] # needs HW mu-maps\n",
+    "except (NameError, KeyError):\n",
+    "    mu = mu_o['im']\n",
+    "    mu_h = {'im': np.zeros_like(mu_o['im'])}"
    ]
   },
   {
@@ -157,10 +167,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "try:  # needs HW mu-maps\n",
-    "    volshow(mu_o['im'] + mu_h['im'], cmaps=['bone'], titles=[r\"$\\mu$-map\"])\n",
-    "except:\n",
-    "    volshow(mu_o['im'], cmaps=['bone'])"
+    "nimpa.imscroll(mu, titles=[r\"$\\mu$-map\"], cmap='bone');"
    ]
   },
   {
@@ -170,11 +177,11 @@
    "outputs": [],
    "source": [
     "# sinogram index (<127 for direct sinograms, >=127 for oblique sinograms)\n",
-    "volshow([m['psino'], m['dsino']],\n",
-    "        titles=[\"Prompt sinogram (%.3gM)\" % (m['psino'].sum() / 1e6),\n",
-    "               \"Delayed sinogram (%.3gM)\" % (m['dsino'].sum() / 1e6)],\n",
-    "        cmaps=['inferno'] * 2, xlabels=[\"\", \"bins\"], ylabels=[\"angles\"] * 2, ncols=2, colorbars=[1]*2,\n",
-    "       figsize=(9.5, 3.5), tight_layout=5);"
+    "nimpa.imscroll(\n",
+    "    {\n",
+    "        f\"Prompt sinogram ({mSino['psino'].sum() / 1e6:.3g}M)\": mSino['psino'],\n",
+    "        f\"Delayed sinogram ({mSino['dsino'].sum() / 1e6:.3g}M)\": mSino['dsino']}, cmap='inferno',\n",
+    "    fig=plt.figure(num=2, figsize=(9.5, 3.5)));"
    ]
   },
   {
@@ -188,7 +195,17 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## OSEM"
+    "## OSEM\n",
+    "\n",
+    "Note that since $\\bf A$ and $\\bf N$ are both diagonal matrix operators, the reconstruction equation can be slightly rewritten to reduce the number of calculations required per iteration:\n",
+    "\n",
+    "$$\n",
+    "{\\bf y}^{(k+1)} = {{\\bf y}^{(k)} \\over \\sum_n{{\\bf H}^T{\\bf X}_n^T{\\bf A}^T{\\bf N}^T{\\bf 1}}}\n",
+    "    \\circ\n",
+    "    \\sum_n{ {\\bf H}^T{\\bf X}_n^T\n",
+    "        { {\\bf m} \\over {\\bf X}_n{\\bf H}{\\bf y}^{(k)} + ({\\bf r} + {\\bf s})/{({\\bf A}^T{\\bf N}^T {\\bf 1})} }\n",
+    "    },\n",
+    "$$"
    ]
   },
   {
@@ -197,20 +214,105 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# built-in default: 14 subsets\n",
-    "recon = nipet.mmrchain(\n",
-    "    datain, mMRpars,\n",
-    "    frames=['timings', [3000, 3600]],\n",
-    "    itr=4,\n",
-    "    histo=m,\n",
-    "    mu_h=mu_h,\n",
-    "    mu_o=mu_o,\n",
-    "    psf='measured',\n",
-    "    outpath=opth,\n",
-    "    fcomment=\"niftypet-recon\",\n",
-    "    store_img=True)\n",
-    "\n",
-    "volshow(recon['im'][:, 100:-100, 100:-100], cmaps=['magma']);"
+    "# setup\n",
+    "psf = nipet.prj.mmrrec.psf_config('measured', mMRpars['Cnt'])\n",
+    "msk = nipet.get_cylinder(mMRpars['Cnt'], rad=29., xo=0., yo=0., unival=1, gpu_dim=True) > 0.9\n",
+    "\n",
+    "## Attenuation, Normalisation & Sensitivity\n",
+    "A = nipet.frwd_prj(mu, mMRpars, attenuation=True, dev_out=True) # imsmooth(mu, psf=psf)\n",
+    "N = nipet.mmrnorm.get_norm_sino(datain, mMRpars, mSino, gpu_dim=True)\n",
+    "AN = cu.asarray(A * N)\n",
+    "\n",
+    "Sn = 14                   # number of subsets\n",
+    "sinoTIdx = [None] * Sn    # subset indices\n",
+    "sen = [None] * Sn         # sensitivity images for each subset\n",
+    "sen_inv_msk = [None] * Sn # masked inverse sensitivity image\n",
+    "\n",
+    "for n in trange(Sn, unit=\"subset\"):\n",
+    "    sinoTIdx[n] = cu.asarray(nipet.prj.mmrrec.get_subsets14(n, mMRpars)[0], 'int32')\n",
+    "    sen[n] = nipet.back_prj(nimpa.isub(AN, sinoTIdx[n]), mMRpars, isub=sinoTIdx[n], dev_out=True,\n",
+    "                            sync=False)\n",
+    "    imsmooth(sen[n], psf=psf, output_array=sen[n])\n",
+    "    assert sen[n].shape == (mMRpars['Cnt']['SZ_IMX'], mMRpars['Cnt']['SZ_IMY'],\n",
+    "                            mMRpars['Cnt']['SZ_IMZ'])\n",
+    "    sen_inv_msk[n] = cu.zeros_like(sen[n])\n",
+    "    sen_inv_msk[n][msk] = np.float32(1) / sen[n][msk]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "## Randoms\n",
+    "\n",
+    "r = nipet.randoms(mSino, mMRpars)[0]\n",
+    "print(\"Randoms: %.3g%%\" % (r.sum() / mSino['psino'].sum() * 100))\n",
+    "\n",
+    "## Scatter\n",
+    "\n",
+    "# Estimated image from two OSEM iterations (implicitly using voxel-driven scatter model)\n",
+    "eim = nipet.mmrchain(datain, mMRpars, mu_h=mu_h, mu_o=mu_o, itr=2, histo=mSino, outpath=opth)['im']\n",
+    "# Recalculate scatter\n",
+    "s = nipet.vsm(datain, (mu_h['im'], mu_o['im']), eim, mMRpars, histo=mSino, rsino=r)\n",
+    "\n",
+    "# normalised randoms + scatter in GPU dimensions\n",
+    "r = nipet.mmraux.remgaps(r, mMRpars['txLUT'], mMRpars['Cnt'])\n",
+    "s = nipet.mmraux.remgaps(s, mMRpars['txLUT'], mMRpars['Cnt'])\n",
+    "rs_AN = nimpa.div(r + s, AN, default=1)\n",
+    "print(\"Scatter: %.3g%%\" % (s.sum() / mSino['psino'].sum() * 100))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "y = cu.ones_like(sen[0])             # reconstructed image\n",
+    "Hy, XHy, crr, bim, mul = (None,) * 5 # temporary variables\n",
+    "m = cu.asarray(nipet.mmraux.remgaps(mSino['psino'], mMRpars['txLUT'], mMRpars['Cnt']))\n",
+    "rs_AN_sub = [nimpa.isub(rs_AN, idx) for idx in sinoTIdx]\n",
+    "\n",
+    "for k in trange(4, desc=\"OSEM\"):\n",
+    "    for n in trange(Sn, unit=\"subset\", leave=False):\n",
+    "        # forward-project current reconstruction to sinogram space\n",
+    "        Hy = imsmooth(y, psf=psf, output_array=Hy, sync=False)\n",
+    "        XHy = nipet.frwd_prj(Hy, mMRpars, isub=sinoTIdx[n], attenuation=False, dev_out=True,\n",
+    "                             fullsino_out=False, output=XHy, sync=False)\n",
+    "        # add randoms and scatter estimates\n",
+    "        if False:\n",
+    "            # recalculate scatter\n",
+    "            s = nipet.vsm(datain, (mu_h['im'], mu_o['im']), y, mMRpars, histo=mSino, rsino=r)\n",
+    "            s = nipet.mmraux.remgaps(s, mMRpars['txLUT'], mMRpars['Cnt'])\n",
+    "            rs_AN = nimpa.div(nimpa.add(nimpa.isub(r, sinoTIdx[n]), nimpa.isub(s, sinoTIdx[n])),\n",
+    "                              nimpa.isub(AN, sinoTIdx[n]), default=1)\n",
+    "            XHy = nimpa.add(XHy, rs_AN, output=XHy, sync=False)\n",
+    "        else:\n",
+    "            XHy = nimpa.add(XHy, rs_AN_sub[n], output=XHy, sync=False)\n",
+    "\n",
+    "        # measured sinogram subset\n",
+    "        crr = nimpa.isub(m, sinoTIdx[n], output=crr, sync=False)\n",
+    "        # corrections\n",
+    "        crr = nimpa.div(crr, XHy, default=1, output=crr, sync=False)\n",
+    "        # back-project corrections to image space\n",
+    "        bim = nipet.back_prj(crr, mMRpars, isub=sinoTIdx[n], dev_out=True, output=bim, sync=False)\n",
+    "        mul = imsmooth(bim, psf=psf, output_array=mul, sync=False)\n",
+    "        # apply FOV mask and normalise with scanner sensitivity\n",
+    "        mul = nimpa.mul(mul, sen_inv_msk[n], output=mul, sync=False)\n",
+    "        # update reconstructed image\n",
+    "        y = nimpa.mul(y, mul, output=y, sync=False)\n",
+    "cu.dev_sync()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nimpa.imscroll(nipet.img.mmrimg.convert2e7(y, mMRpars['Cnt']), cmap='magma', vmin=0, vmax=np.percentile(y, 99.95));"
    ]
   },
   {
@@ -226,25 +328,37 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# setup\n",
+    "psf = nipet.prj.mmrrec.psf_config('measured', mMRpars['Cnt'])\n",
+    "msk = nipet.get_cylinder(mMRpars['Cnt'], rad=29., xo=0., yo=0., unival=1, gpu_dim=True) > 0.9\n",
+    "\n",
     "## Randoms\n",
     "\n",
-    "r = nipet.randoms(m, mMRpars)[0]\n",
-    "print(\"Randoms: %.3g%%\" % (r.sum() / m['psino'].sum() * 100))\n",
+    "r = nipet.randoms(mSino, mMRpars)[0]\n",
+    "print(\"Randoms: %.3g%%\" % (r.sum() / mSino['psino'].sum() * 100))\n",
     "\n",
     "## Scatter\n",
     "\n",
-    "# One OSEM iteration estimate (implicitly using voxel-driven scatter model)\n",
-    "eim = nipet.mmrchain(datain, mMRpars, mu_h=mu_h, mu_o=mu_o, itr=1, histo=m, outpath=opth)['im']\n",
+    "# Estimated image from two OSEM iterations (implicitly using voxel-driven scatter model)\n",
+    "eim = nipet.mmrchain(datain, mMRpars, mu_h=mu_h, mu_o=mu_o, itr=2, histo=mSino, outpath=opth)['im']\n",
     "# Recalculate scatter\n",
-    "s = nipet.vsm(datain, (mu_h['im'], mu_o['im']), eim, mMRpars, histo=m, rsino=r)\n",
-    "print(\"Scatter: %.3g%%\" % (s.sum() / m['psino'].sum() * 100))\n",
+    "s = nipet.vsm(datain, (mu_h['im'], mu_o['im']), eim, mMRpars, histo=mSino, rsino=r)\n",
+    "print(\"Scatter: %.3g%%\" % (s.sum() / mSino['psino'].sum() * 100))\n",
+    "\n",
+    "# normalised randoms + scatter in GPU dimensions\n",
+    "r = nipet.mmraux.remgaps(r, mMRpars['txLUT'], mMRpars['Cnt'])\n",
+    "s = nipet.mmraux.remgaps(s, mMRpars['txLUT'], mMRpars['Cnt'])\n",
     "\n",
     "## Attenuation, Normalisation & Sensitivity\n",
+    "A = nipet.frwd_prj(mu, mMRpars, attenuation=True, dev_out=True)\n",
+    "N = nipet.mmrnorm.get_norm_sino(datain, mMRpars, mSino, gpu_dim=True)\n",
+    "AN = cu.asarray(A * N)\n",
+    "rs_AN = nimpa.div(r + s, AN, default=1)\n",
     "\n",
-    "A = nipet.frwd_prj(mu_h['im'] + mu_o['im'], mMRpars, attenuation=True)\n",
-    "N = nipet.mmrnorm.get_norm_sino(datain, mMRpars, m)\n",
-    "AN = A * N\n",
-    "sim = nipet.back_prj(AN, mMRpars)"
+    "sim = nipet.back_prj(AN, mMRpars, dev_out=True, sync=False)\n",
+    "imsmooth(sim, psf=psf, output_array=sim)\n",
+    "sim_inv_msk = cu.zeros_like(sim)\n",
+    "sim_inv_msk[msk] = np.float32(1) / sim[msk]"
    ]
   },
   {
@@ -253,9 +367,9 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "volshow(OrderedDict([(\"Prompts\", m['psino']), (\"Delayed\", m['dsino']), (\"Attenuation\", A),\n",
-    "                     (\"Scatter\", s), (\"Randoms\", r), (\"Normalisation\", N)]),\n",
-    "        cmaps=['inferno']*6, colorbars=[1]*6, ncols=3, figsize=(9.5, 6));"
+    "# \"Attenuation\": A, \"Normalisation\": N, \"Scatter\": s, \"Randoms\": r\n",
+    "nimpa.imscroll({\"Prompts\": mSino['psino'], \"Delayed\": mSino['dsino']}, cmap='inferno',\n",
+    "               fig=plt.figure(num=4, figsize=(9.5, 3.5)));"
    ]
   },
   {
@@ -265,17 +379,31 @@
    "outputs": [],
    "source": [
     "## MLEM with RM\n",
+    "y = [np.ones_like(sim)]              # reconstructed image\n",
+    "Hy, XHy, crr, bim, mul = (None,) * 5 # temporary variables\n",
+    "m = cu.asarray(nipet.mmraux.remgaps(mSino['psino'], mMRpars['txLUT'], mMRpars['Cnt']))\n",
     "\n",
-    "psf = functools.partial(gaussian_filter, sigma=2.5 / SIGMA2FWHMmm)\n",
-    "msk = nipet.img.mmrimg.get_cylinder(mMRpars['Cnt'], rad=29., xo=0., yo=0., unival=1, gpu_dim=False) <= 0.9\n",
-    "sim_inv = 1 / psf(sim)\n",
-    "sim_inv[msk] = 0\n",
-    "rs = r + s\n",
-    "ANm = AN * m['psino']\n",
-    "recon_mlem = [np.ones_like(sim)]\n",
     "for k in trange(4 * 14, desc=\"MLEM\"):\n",
-    "    fprj = AN * nipet.frwd_prj(psf(recon_mlem[-1]), mMRpars) + rs\n",
-    "    recon_mlem.append(recon_mlem[-1] * sim_inv * psf(nipet.back_prj(ANm / fprj, mMRpars)))"
+    "    # forward-project current reconstruction to sinogram space\n",
+    "    Hy = imsmooth(y[-1], psf=psf, output_array=Hy, sync=False)\n",
+    "    XHy = nipet.frwd_prj(Hy, mMRpars, dev_out=True, fullsino_out=False, output=XHy, sync=False)\n",
+    "    # add randoms and scatter estimates\n",
+    "    if False:\n",
+    "        # recalculate scatter\n",
+    "        s = nipet.vsm(datain, (mu_h['im'], mu_o['im']), y[-1], mMRpars, histo=mSino, rsino=r)\n",
+    "        s = nipet.mmraux.remgaps(s, mMRpars['txLUT'], mMRpars['Cnt'])\n",
+    "        rs_AN = nimpa.div(r + s, AN, default=1)\n",
+    "    XHy = nimpa.add(XHy, rs_AN, output=XHy, sync=False)\n",
+    "\n",
+    "    # corrections\n",
+    "    crr = nimpa.div(m, XHy, default=1, output=crr, sync=False)\n",
+    "    # back-project corrections to image space\n",
+    "    bim = nipet.back_prj(crr, mMRpars, dev_out=True, output=bim, sync=False)\n",
+    "    mul = imsmooth(bim, psf=psf, output_array=mul, sync=False)\n",
+    "    # apply FOV mask and normalise with scanner sensitivity\n",
+    "    mul = nimpa.mul(mul, sim_inv_msk, output=mul, sync=False)\n",
+    "    # update reconstructed image\n",
+    "    y.append(nimpa.mul(y[-1], mul, sync=False))"
    ]
   },
   {
@@ -285,9 +413,11 @@
    "outputs": [],
    "source": [
     "# central slice across iterations\n",
-    "volshow(\n",
-    "    np.asanyarray(recon_mlem[1::5])[:, :, 90:-100, 110:-110],\n",
-    "    cmaps=['magma'] * len(recon_mlem[1::5]), ncols=4, figsize=(9.5, 7), frameon=False);"
+    "nimpa.imscroll(\n",
+    "    {\n",
+    "        f\"{k}\": nipet.img.mmrimg.convert2e7(y[k], mMRpars['Cnt'])[:, 90:-100, 110:-110]\n",
+    "        for k in range(7, len(y), 7)}, cmap='magma', vmin=0, vmax=np.percentile(y[-1], 99.95),\n",
+    "    fig=plt.figure(num=6, figsize=(9.5, 3.5)));"
    ]
   }
  ],
diff --git a/setup.cfg b/setup.cfg
index c538bb46..79333782 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -47,7 +47,7 @@ setup_requires=
     cmake>=3.18
     ninja
 install_requires=
-    cuvec>=2.8.0
+    cuvec>=2.10.0
     miutil>=0.6.0
     nibabel>=2.4.0
     nimpa>=2.0.0

From 0a45341b68bc8b4e512bf2b1365513b46bb1fff2 Mon Sep 17 00:00:00 2001
From: Pawel <p.markiewicz@ucl.ac.uk>
Date: Mon, 3 Jan 2022 17:01:29 +0000
Subject: [PATCH 19/32] added mumap alignment with DIPY and removed pct_mumap
 function

---
 niftypet/nipet/__init__.py   |   2 +-
 niftypet/nipet/img/mmrimg.py | 395 +++++++++++++++++++----------------
 2 files changed, 220 insertions(+), 177 deletions(-)

diff --git a/niftypet/nipet/__init__.py b/niftypet/nipet/__init__.py
index 25c7e116..94ce9ff3 100644
--- a/niftypet/nipet/__init__.py
+++ b/niftypet/nipet/__init__.py
@@ -46,7 +46,7 @@
 from .img.mmrimg import align_mumap
 from .img.mmrimg import convert2dev as im_e72dev
 from .img.mmrimg import convert2e7 as im_dev2e7
-from .img.mmrimg import get_cylinder, hdw_mumap, obj_mumap, pct_mumap
+from .img.mmrimg import get_cylinder, hdw_mumap, obj_mumap
 from .img.pipe import mmrchain
 from .lm.mmrhist import dynamic_timings, mmrhist, randoms
 from .mmraux import explore_input as classify_input
diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index 4b755cc7..8964d06f 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -407,7 +407,7 @@ def align_mumap(
     datain,
     scanner_params=None,
     outpath='',
-    reg_tool='niftyreg',
+    reg_tool='dipy',
     use_stored=False,
     hst=None,
     t0=0,
@@ -572,6 +572,8 @@ def align_mumap(
                 regdct = nimpa.coreg_spm(fpet, fute,
                                          outpath=os.path.join(outpath, 'PET', 'positioning'))
             elif reg_tool == 'niftyreg':
+                if not os.path.exists(Cnt['REGPATH']):
+                    raise ValueError('e> no valid NiftyReg executable')
                 regdct = nimpa.affine_niftyreg(
                     fpet,
                     fute,
@@ -593,6 +595,23 @@ def align_mumap(
                     ffwhm=15.,                                           # pillilitres
                     fthrsh=0.05,
                     verbose=verbose)
+
+            elif reg_tool == 'dipy':
+                regdct = nimpa.affine_dipy(
+                    fpet,
+                    fute,
+                    nbins=32,
+                    metric='MI',
+                    level_iters=[10000, 1000, 200],
+                    sigmas=[3.0, 1.0, 0.0],
+                    factors=[4, 2, 1],
+                    outpath=os.path.join(outpath, 'PET', 'positioning'),
+                    faffine=None,
+                    pickname='ref',
+                    fcomment='',
+                    rfwhm=2.,
+                    ffwhm=2.,
+                    verbose=verbose)
             else:
                 raise ValueError('unknown registration tool requested')
 
@@ -627,6 +646,24 @@ def align_mumap(
                     ffwhm=15.,                                           # pillilitres
                     fthrsh=0.05,
                     verbose=verbose)
+
+            elif reg_tool == 'dipy':
+                regdct = nimpa.affine_dipy(
+                    fpet,
+                    ft1w,
+                    nbins=32,
+                    metric='MI',
+                    level_iters=[10000, 1000, 200],
+                    sigmas=[3.0, 1.0, 0.0],
+                    factors=[4, 2, 1],
+                    outpath=os.path.join(outpath, 'PET', 'positioning'),
+                    faffine=None,
+                    pickname='ref',
+                    fcomment='',
+                    rfwhm=2.,
+                    ffwhm=2.,
+                    verbose=verbose)
+
             else:
                 raise ValueError('unknown registration tool requested')
 
@@ -672,9 +709,11 @@ def align_mumap(
                 raise IOError('The provided NIfTI UTE path is not valid.')
 
     # > call the resampling routine to get the pCT/UTE in place
-    if reg_tool == "spm":
+    if reg_tool == 'spm':
         nimpa.resample_spm(fpet, fflo, faff_mrpet, fimout=freg, del_ref_uncmpr=True,
                            del_flo_uncmpr=True, del_out_uncmpr=True)
+    elif reg_tool == 'dipy':
+        nimpa.resample_dipy(fpet, fflo, faff=faff_mrpet, fimout=freg)
     else:
         nimpa.resample_niftyreg(fpet, fflo, faff_mrpet, fimout=freg, executable=Cnt['RESPATH'],
                                 verbose=verbose)
@@ -732,180 +771,6 @@ def align_mumap(
     return mu_dct
 
 
-# ================================================================================
-# PSEUDO CT MU-MAP
-# --------------------------------------------------------------------------------
-
-
-def pct_mumap(datain, scanner_params, hst=None, t0=0, t1=0, itr=2, petopt='ac', faff='', fpet='',
-              fcomment='', outpath='', store_npy=False, store=False, verbose=False):
-    '''
-    GET THE MU-MAP from pCT IMAGE (which is in T1w space)
-    * the mu-map will be registered to PET which will be reconstructed for time frame t0-t1
-    * it f0 and t1 are not given the whole LM dataset will be reconstructed
-    * the reconstructed PET can be attenuation and scatter corrected or NOT using petopt
-    '''
-    if hst is None:
-        hst = []
-
-    # constants, transaxial and axial LUTs are extracted
-    Cnt = scanner_params['Cnt']
-
-    if not os.path.isfile(faff):
-        from niftypet.nipet.prj import mmrrec
-
-        # histogram the list data if needed
-        if not hst:
-            from niftypet.nipet.lm import mmrhist
-            hst = mmrhist(datain, scanner_params, t0=t0, t1=t1)
-
-    # get hardware mu-map
-    if datain.get("hmumap", "").endswith(".npz") and os.path.isfile(datain["hmumap"]):
-        muh = np.load(datain["hmumap"], allow_pickle=True)["hmu"]
-        (log.info if verbose else log.debug)('loaded hardware mu-map from file:\n{}'.format(
-            datain['hmumap']))
-    elif outpath:
-        hmupath = os.path.join(outpath, "mumap-hdw", "hmumap.npz")
-        if os.path.isfile(hmupath):
-            muh = np.load(hmupath, allow_pickle=True)["hmu"]
-            datain['hmumap'] = hmupath
-        else:
-            raise IOError('Invalid path to the hardware mu-map')
-    else:
-        log.error('The hardware mu-map is required first.')
-        raise IOError('Could not find the hardware mu-map!')
-
-    if not {'MRT1W#', 'T1nii', 'T1bc'}.intersection(datain):
-        log.error('no MR T1w images required for co-registration!')
-        raise IOError('Missing MR data')
-    # ----------------------------------
-
-    # output dictionary
-    mu_dct = {}
-    if not os.path.isfile(faff):
-        # first recon pet to get the T1 aligned to it
-        if petopt == 'qnt':
-            # ---------------------------------------------
-            # OPTION 1 (quantitative recon with all corrections using MR-based mu-map)
-            # get UTE object mu-map (may not be in register with the PET data)
-            mudic = obj_mumap(datain, Cnt)
-            muo = mudic['im']
-            # reconstruct PET image with UTE mu-map to which co-register T1w
-            recout = mmrrec.osemone(datain, [muh, muo], hst, scanner_params, recmod=3, itr=itr,
-                                    fwhm=0., fcomment=fcomment + '_qntUTE',
-                                    outpath=os.path.join(outpath, 'PET',
-                                                         'positioning'), store_img=True)
-        elif petopt == 'nac':
-            # ---------------------------------------------
-            # OPTION 2 (recon without any corrections for scatter and attenuation)
-            # reconstruct PET image with UTE mu-map to which co-register T1w
-            muo = np.zeros(muh.shape, dtype=muh.dtype)
-            recout = mmrrec.osemone(datain, [muh, muo], hst, scanner_params, recmod=1, itr=itr,
-                                    fwhm=0., fcomment=fcomment + '_NAC',
-                                    outpath=os.path.join(outpath, 'PET',
-                                                         'positioning'), store_img=True)
-        elif petopt == 'ac':
-            # ---------------------------------------------
-            # OPTION 3 (recon with attenuation correction only but no scatter)
-            # reconstruct PET image with UTE mu-map to which co-register T1w
-            mudic = obj_mumap(datain, Cnt, outpath=outpath)
-            muo = mudic['im']
-            recout = mmrrec.osemone(datain, [muh, muo], hst, scanner_params, recmod=1, itr=itr,
-                                    fwhm=0., fcomment=fcomment + '_AC',
-                                    outpath=os.path.join(outpath, 'PET',
-                                                         'positioning'), store_img=True)
-
-        fpet = recout.fpet
-        mu_dct['fpet'] = fpet
-
-        # ------------------------------
-        # get the affine transformation
-        ft1w = nimpa.pick_t1w(datain)
-        try:
-            regdct = nimpa.coreg_spm(fpet, ft1w,
-                                     outpath=os.path.join(outpath, 'PET', 'positioning'))
-        except Exception:
-            regdct = nimpa.affine_niftyreg(
-                fpet,
-                ft1w,
-                outpath=os.path.join(outpath, 'PET', 'positioning'), # pcomment=fcomment,
-                executable=Cnt['REGPATH'],
-                omp=multiprocessing.cpu_count() / 2,
-                rigOnly=True,
-                affDirect=False,
-                maxit=5,
-                speed=True,
-                pi=50,
-                pv=50,
-                smof=0,
-                smor=0,
-                rmsk=True,
-                fmsk=True,
-                rfwhm=15.,                                           # pillilitres
-                rthrsh=0.05,
-                ffwhm=15.,                                           # pillilitres
-                fthrsh=0.05,
-                verbose=verbose)
-
-        faff = regdct['faff']
-        # ------------------------------
-
-    # pCT file name
-    if outpath == '':
-        pctdir = os.path.dirname(datain['pCT'])
-    else:
-        pctdir = os.path.join(outpath, 'mumap-obj')
-    mmraux.create_dir(pctdir)
-    fpct = os.path.join(pctdir, 'pCT_r_tmp' + fcomment + '.nii.gz')
-
-    # > call the resampling routine to get the pCT in place
-    if os.path.isfile(Cnt['RESPATH']):
-        cmd = [
-            Cnt['RESPATH'], '-ref', fpet, '-flo', datain['pCT'], '-trans', faff, '-res', fpct,
-            '-pad', '0']
-        if log.getEffectiveLevel() > logging.INFO:
-            cmd.append('-voff')
-        run(cmd)
-    else:
-        log.error('path to resampling executable is incorrect!')
-        raise IOError('Incorrect path to executable!')
-
-    # get the NIfTI of the pCT
-    nim = nib.load(fpct)
-    A = nim.get_sform()
-    pct = nim.get_fdata(dtype=np.float32)
-    pct = pct[:, ::-1, ::-1]
-    pct = np.transpose(pct, (2, 1, 0))
-    # convert the HU units to mu-values
-    mu = hu2mu(pct)
-    # get rid of negatives
-    mu[mu < 0] = 0
-
-    # return image dictionary with the image itself and other parameters
-    mu_dct['im'] = mu
-    mu_dct['affine'] = A
-    mu_dct['faff'] = faff
-
-    if store:
-        # now save to numpy array and NIfTI in this folder
-        if outpath == '':
-            pctumapdir = os.path.join(datain['corepath'], 'mumap-obj')
-        else:
-            pctumapdir = os.path.join(outpath, 'mumap-obj')
-        mmraux.create_dir(pctumapdir)
-        # > Numpy
-        if store_npy:
-            fnp = os.path.join(pctumapdir, "mumap-pCT.npz")
-            np.savez(fnp, mu=mu, A=A)
-
-        # > NIfTI
-        fmu = os.path.join(pctumapdir, 'mumap-pCT' + fcomment + '.nii.gz')
-        nimpa.array2nii(mu[::-1, ::-1, :], A, fmu)
-        mu_dct['fim'] = fmu
-        datain['mumapCT'] = fmu
-
-    return mu_dct
-
 
 # ********************************************************************************
 # GET HARDWARE MU-MAPS with positions and offsets
@@ -1314,3 +1179,181 @@ def rmumaps(datain, Cnt, t0=0, t1=0, use_stored=False):
         muh = muh[2 * Cnt['RNG_STRT']:2 * Cnt['RNG_END'], :, :]
         muo = muo[2 * Cnt['RNG_STRT']:2 * Cnt['RNG_END'], :, :]
         return [muh, muo]
+
+
+
+
+
+# # ================================================================================
+# # PSEUDO CT MU-MAP
+# # --------------------------------------------------------------------------------
+
+
+# def pct_mumap(datain, scanner_params, hst=None, t0=0, t1=0, itr=2, petopt='ac', faff='', fpet='',
+#               fcomment='', outpath='', store_npy=False, store=False, verbose=False):
+#     '''
+#     GET THE MU-MAP from pCT IMAGE (which is in T1w space)
+#     * the mu-map will be registered to PET which will be reconstructed for time frame t0-t1
+#     * it f0 and t1 are not given the whole LM dataset will be reconstructed
+#     * the reconstructed PET can be attenuation and scatter corrected or NOT using petopt
+#     '''
+#     if hst is None:
+#         hst = []
+
+#     # constants, transaxial and axial LUTs are extracted
+#     Cnt = scanner_params['Cnt']
+
+#     if not os.path.isfile(faff):
+#         from niftypet.nipet.prj import mmrrec
+
+#         # histogram the list data if needed
+#         if not hst:
+#             from niftypet.nipet.lm import mmrhist
+#             hst = mmrhist(datain, scanner_params, t0=t0, t1=t1)
+
+#     # get hardware mu-map
+#     if datain.get("hmumap", "").endswith(".npz") and os.path.isfile(datain["hmumap"]):
+#         muh = np.load(datain["hmumap"], allow_pickle=True)["hmu"]
+#         (log.info if verbose else log.debug)('loaded hardware mu-map from file:\n{}'.format(
+#             datain['hmumap']))
+#     elif outpath:
+#         hmupath = os.path.join(outpath, "mumap-hdw", "hmumap.npz")
+#         if os.path.isfile(hmupath):
+#             muh = np.load(hmupath, allow_pickle=True)["hmu"]
+#             datain['hmumap'] = hmupath
+#         else:
+#             raise IOError('Invalid path to the hardware mu-map')
+#     else:
+#         log.error('The hardware mu-map is required first.')
+#         raise IOError('Could not find the hardware mu-map!')
+
+#     if not {'MRT1W#', 'T1nii', 'T1bc'}.intersection(datain):
+#         log.error('no MR T1w images required for co-registration!')
+#         raise IOError('Missing MR data')
+#     # ----------------------------------
+
+#     # output dictionary
+#     mu_dct = {}
+#     if not os.path.isfile(faff):
+#         # first recon pet to get the T1 aligned to it
+#         if petopt == 'qnt':
+#             # ---------------------------------------------
+#             # OPTION 1 (quantitative recon with all corrections using MR-based mu-map)
+#             # get UTE object mu-map (may not be in register with the PET data)
+#             mudic = obj_mumap(datain, Cnt)
+#             muo = mudic['im']
+#             # reconstruct PET image with UTE mu-map to which co-register T1w
+#             recout = mmrrec.osemone(datain, [muh, muo], hst, scanner_params, recmod=3, itr=itr,
+#                                     fwhm=0., fcomment=fcomment + '_qntUTE',
+#                                     outpath=os.path.join(outpath, 'PET',
+#                                                          'positioning'), store_img=True)
+#         elif petopt == 'nac':
+#             # ---------------------------------------------
+#             # OPTION 2 (recon without any corrections for scatter and attenuation)
+#             # reconstruct PET image with UTE mu-map to which co-register T1w
+#             muo = np.zeros(muh.shape, dtype=muh.dtype)
+#             recout = mmrrec.osemone(datain, [muh, muo], hst, scanner_params, recmod=1, itr=itr,
+#                                     fwhm=0., fcomment=fcomment + '_NAC',
+#                                     outpath=os.path.join(outpath, 'PET',
+#                                                          'positioning'), store_img=True)
+#         elif petopt == 'ac':
+#             # ---------------------------------------------
+#             # OPTION 3 (recon with attenuation correction only but no scatter)
+#             # reconstruct PET image with UTE mu-map to which co-register T1w
+#             mudic = obj_mumap(datain, Cnt, outpath=outpath)
+#             muo = mudic['im']
+#             recout = mmrrec.osemone(datain, [muh, muo], hst, scanner_params, recmod=1, itr=itr,
+#                                     fwhm=0., fcomment=fcomment + '_AC',
+#                                     outpath=os.path.join(outpath, 'PET',
+#                                                          'positioning'), store_img=True)
+
+#         fpet = recout.fpet
+#         mu_dct['fpet'] = fpet
+
+#         # ------------------------------
+#         # get the affine transformation
+#         ft1w = nimpa.pick_t1w(datain)
+#         try:
+#             regdct = nimpa.coreg_spm(fpet, ft1w,
+#                                      outpath=os.path.join(outpath, 'PET', 'positioning'))
+#         except Exception:
+#             regdct = nimpa.affine_niftyreg(
+#                 fpet,
+#                 ft1w,
+#                 outpath=os.path.join(outpath, 'PET', 'positioning'), # pcomment=fcomment,
+#                 executable=Cnt['REGPATH'],
+#                 omp=multiprocessing.cpu_count() / 2,
+#                 rigOnly=True,
+#                 affDirect=False,
+#                 maxit=5,
+#                 speed=True,
+#                 pi=50,
+#                 pv=50,
+#                 smof=0,
+#                 smor=0,
+#                 rmsk=True,
+#                 fmsk=True,
+#                 rfwhm=15.,                                           # pillilitres
+#                 rthrsh=0.05,
+#                 ffwhm=15.,                                           # pillilitres
+#                 fthrsh=0.05,
+#                 verbose=verbose)
+
+#         faff = regdct['faff']
+#         # ------------------------------
+
+#     # pCT file name
+#     if outpath == '':
+#         pctdir = os.path.dirname(datain['pCT'])
+#     else:
+#         pctdir = os.path.join(outpath, 'mumap-obj')
+#     mmraux.create_dir(pctdir)
+#     fpct = os.path.join(pctdir, 'pCT_r_tmp' + fcomment + '.nii.gz')
+
+#     # > call the resampling routine to get the pCT in place
+#     if os.path.isfile(Cnt['RESPATH']):
+#         cmd = [
+#             Cnt['RESPATH'], '-ref', fpet, '-flo', datain['pCT'], '-trans', faff, '-res', fpct,
+#             '-pad', '0']
+#         if log.getEffectiveLevel() > logging.INFO:
+#             cmd.append('-voff')
+#         run(cmd)
+#     else:
+#         log.error('path to resampling executable is incorrect!')
+#         raise IOError('Incorrect path to executable!')
+
+#     # get the NIfTI of the pCT
+#     nim = nib.load(fpct)
+#     A = nim.get_sform()
+#     pct = nim.get_fdata(dtype=np.float32)
+#     pct = pct[:, ::-1, ::-1]
+#     pct = np.transpose(pct, (2, 1, 0))
+#     # convert the HU units to mu-values
+#     mu = hu2mu(pct)
+#     # get rid of negatives
+#     mu[mu < 0] = 0
+
+#     # return image dictionary with the image itself and other parameters
+#     mu_dct['im'] = mu
+#     mu_dct['affine'] = A
+#     mu_dct['faff'] = faff
+
+#     if store:
+#         # now save to numpy array and NIfTI in this folder
+#         if outpath == '':
+#             pctumapdir = os.path.join(datain['corepath'], 'mumap-obj')
+#         else:
+#             pctumapdir = os.path.join(outpath, 'mumap-obj')
+#         mmraux.create_dir(pctumapdir)
+#         # > Numpy
+#         if store_npy:
+#             fnp = os.path.join(pctumapdir, "mumap-pCT.npz")
+#             np.savez(fnp, mu=mu, A=A)
+
+#         # > NIfTI
+#         fmu = os.path.join(pctumapdir, 'mumap-pCT' + fcomment + '.nii.gz')
+#         nimpa.array2nii(mu[::-1, ::-1, :], A, fmu)
+#         mu_dct['fim'] = fmu
+#         datain['mumapCT'] = fmu
+
+#     return mu_dct

From 11047a26e724de48d3592895d0407f8638e00fb9 Mon Sep 17 00:00:00 2001
From: Pawel <p.markiewicz@ucl.ac.uk>
Date: Mon, 3 Jan 2022 21:12:40 +0000
Subject: [PATCH 20/32] small fixes in mumap alignment

---
 niftypet/nipet/img/mmrimg.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index 8964d06f..dc6f8b2d 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -625,6 +625,8 @@ def align_mumap(
                 regdct = nimpa.coreg_spm(fpet, ft1w,
                                          outpath=os.path.join(outpath, 'PET', 'positioning'))
             elif reg_tool == 'niftyreg':
+                if not os.path.exists(Cnt['REGPATH']):
+                    raise ValueError('e> no valid NiftyReg executable')
                 regdct = nimpa.affine_niftyreg(
                     fpet,
                     ft1w,

From c69e0cab83e724ebd3c92ba46f4629d1fe3d3fc7 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sun, 2 Jan 2022 11:48:15 +0000
Subject: [PATCH 21/32] lint

---
 niftypet/nipet/img/mmrimg.py | 28 +++++++---------------------
 niftypet/nipet/img/pipe.py   | 13 +++----------
 niftypet/nipet/prj/mmrrec.py |  3 ++-
 3 files changed, 12 insertions(+), 32 deletions(-)

diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index dc6f8b2d..c363be73 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -278,15 +278,10 @@ def mudcm2nii(datain, Cnt):
     im = np.zeros((Cnt['SO_IMZ'], Cnt['SO_IMY'], Cnt['SO_IMX']), dtype=np.float32)
     nimpa.array2nii(im, B, os.path.join(os.path.dirname(datain['mumapDCM']), 'muref.nii.gz'))
     # -------------------------------------------------------------------------------------
-    fmu = os.path.join(os.path.dirname(datain['mumapDCM']), 'mu_r.nii.gz')
-
-    nimpa.resample_dipy(
-        os.path.join(os.path.dirname(datain['mumapDCM']), 'muref.nii.gz'),
-        os.path.join(os.path.dirname(datain['mumapDCM']), 'mu.nii.gz'),
-        fimout=fmu,
-        intrp=1,
-        dtype_nifti=np.float32)
-
+    opth = os.path.dirname(datain['mumapDCM'])
+    fmu = os.path.join(opth, 'mu_r.nii.gz')
+    nimpa.resample_dipy(os.path.join(opth, 'muref.nii.gz'), os.path.join(opth, 'mu.nii.gz'),
+                        fimout=fmu, intrp=1, dtype_nifti=np.float32)
     return fmu
 
 
@@ -353,12 +348,7 @@ def obj_mumap(
 
     # > resampled the NIfTI converted image to the reference shape/size
     fmu = os.path.join(fmudir, comment + 'mumap_tmp.nii.gz')
-    nimpa.resample_dipy(
-        fmuref,
-        fmunii,
-        fimout=fmu,
-        intrp=1,
-        dtype_nifti=np.float32)
+    nimpa.resample_dipy(fmuref, fmunii, fimout=fmu, intrp=1, dtype_nifti=np.float32)
 
     nim = nib.load(fmu)
     # get the affine transform
@@ -1007,12 +997,8 @@ def get_hmupos(datain, parts, Cnt, outpath=''):
         fout = os.path.join(os.path.dirname(hmupos[0]['niipath']),
                             'r' + os.path.basename(hmupos[i]['niipath']).split('.')[0] + '.nii.gz')
 
-        nimpa.resample_dipy(
-            hmupos[0]['niipath'],
-            hmupos[i]['niipath'],
-            fimout=fout,
-            intrp=1,
-            dtype_nifti=np.float32)
+        nimpa.resample_dipy(hmupos[0]['niipath'], hmupos[i]['niipath'], fimout=fout, intrp=1,
+                            dtype_nifti=np.float32)
 
     return hmupos
 
diff --git a/niftypet/nipet/img/pipe.py b/niftypet/nipet/img/pipe.py
index 5ca1b8dd..bff45dc4 100644
--- a/niftypet/nipet/img/pipe.py
+++ b/niftypet/nipet/img/pipe.py
@@ -48,7 +48,6 @@ def mmrchain(
     trim_memlim=True,       # reduced use of memory for machines
                             # with limited memory (slow though)
     trim_store=True,
-
     pvcroi=None,            # ROI used for PVC.  If undefined no PVC
                             # is performed.
     pvcreg_tool='niftyreg', # the registration tool used in PVC
@@ -408,15 +407,9 @@ def mmrchain(
         elif 'lm_ima' in datain:
             fnm = os.path.basename(datain['lm_ima'])[:20]
         # trim PET and upsample
-        petu = nimpa.imtrimup(dynim,
-                              affine=image_affine(datain, Cnt),
-                              scale=trim_scale,
-                              int_order=trim_interp,
-                              outpath=petimg,
-                              fname=fnm,
-                              fcomment=fcomment,
-                              fcomment_pfx=fout+'_',
-                              store_img=trim_store, 
+        petu = nimpa.imtrimup(dynim, affine=image_affine(datain, Cnt), scale=trim_scale,
+                              int_order=trim_interp, outpath=petimg, fname=fnm, fcomment=fcomment,
+                              fcomment_pfx=fout + '_', store_img=trim_store,
                               store_img_intrmd=store_img_intrmd, memlim=trim_memlim,
                               verbose=log.getEffectiveLevel())
 
diff --git a/niftypet/nipet/prj/mmrrec.py b/niftypet/nipet/prj/mmrrec.py
index 841469cd..a2b2411d 100644
--- a/niftypet/nipet/prj/mmrrec.py
+++ b/niftypet/nipet/prj/mmrrec.py
@@ -232,7 +232,8 @@ def osemone(datain, mumaps, hst, scanner_params, recmod=3, itr=4, fwhm=0., psf=N
             log.info('using provided attenuation factor sinogram')
         else:
             asng = cu.zeros(psng.shape, dtype=np.float32)
-            petprj.fprj(asng, cu.asarray(mus), txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt, 1)
+            petprj.fprj(asng, cu.asarray(mus), txLUT, axLUT, np.array([-1], dtype=np.int32), Cnt,
+                        1)
     # > combine attenuation and normalisation
     ansng = asng * nsng
     # ========================================================================

From d8d6b981a0fa67d4f71b033a428efc54da74c371 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Tue, 8 Feb 2022 02:19:00 +0000
Subject: [PATCH 22/32] use nimpa.dcm2nii

---
 niftypet/nipet/img/mmrimg.py | 68 ++++++++++--------------------------
 1 file changed, 18 insertions(+), 50 deletions(-)

diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index c363be73..d5f07e41 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -340,11 +340,7 @@ def obj_mumap(
         os.remove(d)
 
     # > convert the DICOM mu-map images to NIfTI
-    run([Cnt['DCM2NIIX'], '-f', fnii + tstmp, '-o', fmudir, datain['mumapDCM']])
-    # piles for the T1w, pick one:
-    fmunii = glob.glob(os.path.join(fmudir, '*' + fnii + tstmp + '*.nii*'))[0]
-    # fmunii = glob.glob( os.path.join(datain['mumapDCM'], '*converted*.nii*') )
-    # fmunii = fmunii[0]
+    fmunii = nimpa.dcm2nii(datain['mumapDCM'], fnii + tstmp, outpath=fmudir)
 
     # > resampled the NIfTI converted image to the reference shape/size
     fmu = os.path.join(fmudir, comment + 'mumap_tmp.nii.gz')
@@ -551,9 +547,7 @@ def align_mumap(
         if musrc == 'ute' and ute_name in datain and os.path.exists(datain[ute_name]):
             # change to NIfTI if the UTE sequence is in DICOM files (folder)
             if os.path.isdir(datain[ute_name]):
-                fnew = os.path.basename(datain[ute_name])
-                run([Cnt['DCM2NIIX'], '-f', fnew, datain[ute_name]])
-                fute = glob.glob(os.path.join(datain[ute_name], fnew + '*nii*'))[0]
+                fute = nimpa.dcm2nii(datain[ute_name], os.path.basename(datain[ute_name]))
             elif os.path.isfile(datain[ute_name]):
                 fute = datain[ute_name]
 
@@ -587,21 +581,13 @@ def align_mumap(
                     verbose=verbose)
 
             elif reg_tool == 'dipy':
-                regdct = nimpa.affine_dipy(
-                    fpet,
-                    fute,
-                    nbins=32,
-                    metric='MI',
-                    level_iters=[10000, 1000, 200],
-                    sigmas=[3.0, 1.0, 0.0],
-                    factors=[4, 2, 1],
-                    outpath=os.path.join(outpath, 'PET', 'positioning'),
-                    faffine=None,
-                    pickname='ref',
-                    fcomment='',
-                    rfwhm=2.,
-                    ffwhm=2.,
-                    verbose=verbose)
+                regdct = nimpa.affine_dipy(fpet, fute, nbins=32, metric='MI',
+                                           level_iters=[10000, 1000,
+                                                        200], sigmas=[3.0, 1.0,
+                                                                      0.0], factors=[4, 2, 1],
+                                           outpath=os.path.join(outpath, 'PET', 'positioning'),
+                                           faffine=None, pickname='ref', fcomment='', rfwhm=2.,
+                                           ffwhm=2., verbose=verbose)
             else:
                 raise ValueError('unknown registration tool requested')
 
@@ -640,21 +626,13 @@ def align_mumap(
                     verbose=verbose)
 
             elif reg_tool == 'dipy':
-                regdct = nimpa.affine_dipy(
-                    fpet,
-                    ft1w,
-                    nbins=32,
-                    metric='MI',
-                    level_iters=[10000, 1000, 200],
-                    sigmas=[3.0, 1.0, 0.0],
-                    factors=[4, 2, 1],
-                    outpath=os.path.join(outpath, 'PET', 'positioning'),
-                    faffine=None,
-                    pickname='ref',
-                    fcomment='',
-                    rfwhm=2.,
-                    ffwhm=2.,
-                    verbose=verbose)
+                regdct = nimpa.affine_dipy(fpet, ft1w, nbins=32, metric='MI',
+                                           level_iters=[10000, 1000,
+                                                        200], sigmas=[3.0, 1.0,
+                                                                      0.0], factors=[4, 2, 1],
+                                           outpath=os.path.join(outpath, 'PET', 'positioning'),
+                                           faffine=None, pickname='ref', fcomment='', rfwhm=2.,
+                                           ffwhm=2., verbose=verbose)
 
             else:
                 raise ValueError('unknown registration tool requested')
@@ -691,9 +669,7 @@ def align_mumap(
             # convert the DICOM mu-map images to nii
             if 'mumapDCM' not in datain:
                 raise IOError('DICOM with the UTE mu-map are not given.')
-            run([Cnt['DCM2NIIX'], '-f', fnii + tstmp, '-o', opth, datain['mumapDCM']])
-            # piles for the T1w, pick one:
-            fflo = glob.glob(os.path.join(opth, '*' + fnii + tstmp + '*.nii*'))[0]
+            fflo = nimpa.dcm2nii(datain['mumapDCM'], fnii + tstmp, outpath=opth)
         else:
             if os.path.isfile(datain['UTE']):
                 fflo = datain['UTE']
@@ -763,7 +739,6 @@ def align_mumap(
     return mu_dct
 
 
-
 # ********************************************************************************
 # GET HARDWARE MU-MAPS with positions and offsets
 # --------------------------------------------------------------------------------
@@ -1132,10 +1107,7 @@ def rmumaps(datain, Cnt, t0=0, t1=0, use_stored=False):
             ft1w = datain['T1bc']
         elif os.path.isdir(datain['MRT1W']):
             # create file name for the converted NIfTI image
-            fnii = 'converted'
-            run([Cnt['DCM2NIIX'], '-f', fnii, datain['T1nii']])
-            ft1nii = glob.glob(os.path.join(datain['T1nii'], '*converted*.nii*'))
-            ft1w = ft1nii[0]
+            ft1w = nimpa.dcm2nii(datain['T1nii'], 'converted')
         else:
             raise IOError('Disaster: no T1w image!')
 
@@ -1169,14 +1141,10 @@ def rmumaps(datain, Cnt, t0=0, t1=0, use_stored=False):
         return [muh, muo]
 
 
-
-
-
 # # ================================================================================
 # # PSEUDO CT MU-MAP
 # # --------------------------------------------------------------------------------
 
-
 # def pct_mumap(datain, scanner_params, hst=None, t0=0, t1=0, itr=2, petopt='ac', faff='', fpet='',
 #               fcomment='', outpath='', store_npy=False, store=False, verbose=False):
 #     '''

From 311ad0f4be712e6aeea193ea82db5de00bae0c19 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Tue, 8 Feb 2022 19:43:45 +0000
Subject: [PATCH 23/32] img/pipe: fix dynamic_timings

- fixes #43
---
 niftypet/nipet/img/pipe.py | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/niftypet/nipet/img/pipe.py b/niftypet/nipet/img/pipe.py
index bff45dc4..89df9b31 100644
--- a/niftypet/nipet/img/pipe.py
+++ b/niftypet/nipet/img/pipe.py
@@ -113,15 +113,11 @@ def mmrchain(
         elif (isinstance(frames[0], str) and frames[0] == 'def'
               and all(isinstance(t, list) and len(t) == 2 for t in frames[1:])):
             # get total time and list of all time frames
-            dfrms = dynamic_timings(frames)
-            t_frms = dfrms[1:]
-
+            t_frms = dynamic_timings(frames)['timings'][1:]
         # if 1D:
         elif all(isinstance(t, Integral) for t in frames):
             # get total time and list of all time frames
-            dfrms = dynamic_timings(frames)
-            t_frms = dfrms[1:]
-
+            t_frms = dynamic_timings(frames)['timings'][1:]
         else:
             log.error('osemdyn: frames definitions are not given\
                 in the correct list format: 1D [15,15,30,30,...]\

From 11cf97ca600177214cf26332dff4b10c18355eb6 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Tue, 8 Feb 2022 19:44:55 +0000
Subject: [PATCH 24/32] dynamic_timings: sync img/auximg with lm/mmrhist

---
 niftypet/nipet/img/auximg.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/niftypet/nipet/img/auximg.py b/niftypet/nipet/img/auximg.py
index 3cfbce4b..9e59f14b 100644
--- a/niftypet/nipet/img/auximg.py
+++ b/niftypet/nipet/img/auximg.py
@@ -76,9 +76,10 @@ def dynamic_timings(flist, offset=0):
     Arguments:
       flist: can be 1D list of time duration for each dynamic frame, e.g.:
             flist = [15, 15, 15, 15, 30, 30, 30, ...]
-        or a 2D list of lists having 2 entries:
+        or a 2D list of lists having 2 entries per definition:
         first for the number of repetitions and the other for the frame duration, e.g.:
-            flist = [[4,15], [3,15], ...].
+            flist = [[4, 15], [8, 30], ...],
+        meaning 4x15s, then 8x30s, etc.
       offset: adjusts for the start time (usually when prompts are strong enough over randoms)
     Returns (dict):
       'timings': [[0, 15], [15, 30], [30, 45], [45, 60], [60, 90], [90, 120], [120, 150], ...]
@@ -119,8 +120,8 @@ def dynamic_timings(flist, offset=0):
         tsum = 0
         # list of frame timings
         t_frames = []
-        for i in range(0, farray.shape[0]):
-            for _ in range(0, farray[i, 0]):
+        for i in range(farray.shape[0]):
+            for _ in range(farray[i, 0]):
                 # frame start time
                 t0 = tsum
                 tsum += farray[i, 1]
@@ -131,5 +132,5 @@ def dynamic_timings(flist, offset=0):
                 frms[fi] = farray[i, 1]
                 fi += 1
     else:
-        raise TypeError('Unrecognised data input.')
+        raise TypeError('Unrecognised time frame definitions.')
     return {'total': tsum, 'frames': frms, 'timings': t_frames}

From 24661214d43fbfc4f7938ac55bd8f10ced48a7eb Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Fri, 20 May 2022 08:36:14 +0100
Subject: [PATCH 25/32] linting & pre-commit updates

---
 .pre-commit-config.yaml      |  7 +++++--
 niftypet/nipet/img/mmrimg.py | 14 ++++++--------
 niftypet/nipet/img/pipe.py   |  5 ++---
 niftypet/nipet/mmraux.py     |  4 ++--
 niftypet/nipet/prj/mmrrec.py | 12 ++++++------
 niftypet/nipet/prj/mmrsim.py |  6 ++----
 6 files changed, 23 insertions(+), 25 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 77e4b347..929980db 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.2.0
   hooks:
   - id: check-added-large-files
   - id: check-case-conflict
@@ -31,15 +31,18 @@ repos:
   - id: flake8
     args: [-j8]
     additional_dependencies:
+    - flake8-broken-line
     - flake8-bugbear
     - flake8-comprehensions
     - flake8-debugger
+    - 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:
diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index d5f07e41..c201f8a1 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -442,9 +442,8 @@ def align_mumap(
     # ---------------------------------------------------------------------------
     # > used stored if requested
     if use_stored:
-        fmu_stored = fnm + '-aligned-to_t'\
-                     + str(t0)+'-'+str(t1)+'_'+petopt.upper()\
-                     + fcomment
+        fmu_stored = fnm + '-aligned-to_t' + str(t0) + '-' + str(
+            t1) + '_' + petopt.upper() + fcomment
         fmupath = os.path.join(opth, fmu_stored + '.nii.gz')
 
         if os.path.isfile(fmupath):
@@ -478,8 +477,8 @@ def align_mumap(
             if 'txLUT' in scanner_params:
                 hst = mmrhist(datain, scanner_params, t0=t0, t1=t1)
             else:
-                raise ValueError('Full scanner are parameters not provided\
-                     but are required for histogramming.')
+                raise ValueError(
+                    'Full scanner are parameters not provided but are required for histogramming.')
 
     # ========================================================
     # -get hardware mu-map
@@ -715,9 +714,8 @@ def align_mumap(
     if store or store_npy:
         nimpa.create_dir(opth)
         if faff == '':
-            fname = fnm + '-aligned-to_t'\
-                    + str(t0)+'-'+str(t1)+'_'+petopt.upper()\
-                    + fcomment
+            fname = fnm + '-aligned-to_t' + str(t0) + '-' + str(
+                t1) + '_' + petopt.upper() + fcomment
         else:
             fname = fnm + '-aligned-to-given-affine' + fcomment
     if store_npy:
diff --git a/niftypet/nipet/img/pipe.py b/niftypet/nipet/img/pipe.py
index 89df9b31..351d4206 100644
--- a/niftypet/nipet/img/pipe.py
+++ b/niftypet/nipet/img/pipe.py
@@ -119,9 +119,8 @@ def mmrchain(
             # get total time and list of all time frames
             t_frms = dynamic_timings(frames)['timings'][1:]
         else:
-            log.error('osemdyn: frames definitions are not given\
-                in the correct list format: 1D [15,15,30,30,...]\
-                or 2D list [[2,15], [2,30], ...]')
+            log.error('osemdyn: frames definitions are not given in the correct list format:'
+                      ' 1D [15,15,30,30,...] or 2D list [[2,15], [2,30], ...]')
     else:
         log.error(
             'provided dynamic frames definitions are incorrect (should be a list of definitions).')
diff --git a/niftypet/nipet/mmraux.py b/niftypet/nipet/mmraux.py
index a70f4f13..54584220 100644
--- a/niftypet/nipet/mmraux.py
+++ b/niftypet/nipet/mmraux.py
@@ -1023,7 +1023,7 @@ def get_dicoms(dfile, datain, Cnt):
         else:
             f0 = csahdr.find('RadionuclideCodeSequence')
             if f0 < 0:
-                print('w> could not find isotope name.  enter manually into Cnt[' 'ISOTOPE' ']')
+                log.warning("Could not find isotope name. Enter manually into Cnt['ISOTOPE']")
                 return None
             istp_coded = re.search(r'(?<=CodeValue:)\S*', csahdr[f0:f0 + 100]).group()
             if istp_coded == 'C-111A1': Cnt['ISOTOPE'] = 'F18'
@@ -1032,7 +1032,7 @@ def get_dicoms(dfile, datain, Cnt):
             elif istp_coded == 'C-128A2': Cnt['ISOTOPE'] = 'Ge68'
             elif istp_coded == 'C-131A3': Cnt['ISOTOPE'] = 'Ga68'
             else:
-                print('w> could not find isotope name.  enter manually into Cnt[' 'ISOTOPE' ']')
+                log.warning("Could not find isotope name. Enter manually into Cnt['ISOTOPE']")
                 return None
         # ---
 
diff --git a/niftypet/nipet/prj/mmrrec.py b/niftypet/nipet/prj/mmrrec.py
index a2b2411d..62cd7f99 100644
--- a/niftypet/nipet/prj/mmrrec.py
+++ b/niftypet/nipet/prj/mmrrec.py
@@ -105,6 +105,7 @@ def psf_config(psf, Cnt):
         [x, y, z]: list or Numpy array of separate FWHM of the PSF for each direction
         ndarray: 3 x 2*RSZ_PSF_KRNL+1 Numpy array directly defining the kernel in each direction
     '''
+
     def _config(fwhm3, check_len=True):
         # resolution modelling by custom kernels
         if check_len:
@@ -222,12 +223,11 @@ def osemone(datain, mumaps, hst, scanner_params, recmod=3, itr=4, fwhm=0., psf=N
         asng = np.ones(psng.shape, dtype=np.float32)
     else:
         # > check if the attenuation sino is given as an array
-        if isinstance(attnsino, np.ndarray) \
-                and attnsino.shape==(Cnt['NSN11'], Cnt['NSANGLES'], Cnt['NSBINS']):
+        if isinstance(attnsino, np.ndarray) and attnsino.shape == (Cnt['NSN11'], Cnt['NSANGLES'],
+                                                                   Cnt['NSBINS']):
             asng = mmraux.remgaps(attnsino, txLUT, Cnt)
             log.info('using provided attenuation factor sinogram')
-        elif isinstance(attnsino, np.ndarray) \
-                and attnsino.shape==(Cnt['Naw'], Cnt['NSN11']):
+        elif isinstance(attnsino, np.ndarray) and attnsino.shape == (Cnt['Naw'], Cnt['NSN11']):
             asng = attnsino
             log.info('using provided attenuation factor sinogram')
         else:
@@ -241,8 +241,8 @@ def osemone(datain, mumaps, hst, scanner_params, recmod=3, itr=4, fwhm=0., psf=N
     # ========================================================================
     # Randoms
     # -------------------------------------------------------------------------
-    if isinstance(randsino, np.ndarray) \
-            and randsino.shape==(Cnt['NSN11'], Cnt['NSANGLES'], Cnt['NSBINS']):
+    if isinstance(randsino,
+                  np.ndarray) and randsino.shape == (Cnt['NSN11'], Cnt['NSANGLES'], Cnt['NSBINS']):
         rsino = randsino
         rsng = mmraux.remgaps(randsino, txLUT, Cnt)
     else:
diff --git a/niftypet/nipet/prj/mmrsim.py b/niftypet/nipet/prj/mmrsim.py
index 1ae5eb3f..f84c0208 100644
--- a/niftypet/nipet/prj/mmrsim.py
+++ b/niftypet/nipet/prj/mmrsim.py
@@ -43,8 +43,7 @@ def simulate_sino(
         raise ValueError('The shapes of the PET and CT images are inconsistent.')
 
     if simulate_3d:
-        if petim.ndim != 3 \
-                or petim.shape != (Cnt['SO_IMZ'], Cnt['SO_IMY'], Cnt['SO_IMX']):
+        if petim.ndim != 3 or petim.shape != (Cnt['SO_IMZ'], Cnt['SO_IMY'], Cnt['SO_IMX']):
             raise ValueError('The input image shape does not match the scanner image size.')
         if petim.max() > 200:
             log.warning('the PET image may have too large intensities for robust simulation.')
@@ -152,8 +151,7 @@ def simulate_recon(
     psfkernel = mmrrec.psf_config(psf, Cnt)
 
     if simulate_3d:
-        if ctim.ndim!=3 \
-                or ctim.shape!=(Cnt['SO_IMZ'], Cnt['SO_IMY'], Cnt['SO_IMX']):
+        if ctim.ndim != 3 or ctim.shape != (Cnt['SO_IMZ'], Cnt['SO_IMY'], Cnt['SO_IMX']):
             raise ValueError('The CT/mu-map image does not match the scanner image shape.')
     else:
         # > 2D case with reduced rings

From 811666066656fb64620acaf92f399f94b311e4b8 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Fri, 20 May 2022 09:00:50 +0100
Subject: [PATCH 26/32] bump nimpa, drop regpath/respath

---
 .github/workflows/test.yml   |  3 +--
 niftypet/nipet/img/mmrimg.py |  9 +--------
 niftypet/nipet/img/pipe.py   | 13 ++-----------
 setup.cfg                    |  2 +-
 4 files changed, 5 insertions(+), 22 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 2593fa48..7aee93e2 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -39,7 +39,7 @@ jobs:
     runs-on: [self-hosted, python, cuda, matlab]
     strategy:
       matrix:
-        python: [3.7, 3.8]
+        python: [3.7, 3.9]
     steps:
     - uses: actions/checkout@v2
       with:
@@ -70,7 +70,6 @@ jobs:
         password: ${{ secrets.PYPI_TOKEN }}
         upload: ${{ github.event_name == 'push' && startsWith(github.ref, 'refs/tags') }}
       env:
-        PATHTOOLS: ${{ github.workspace }}/NiftyPET_tools
         HMUDIR: ${{ github.workspace }}
     - id: meta
       name: Changelog
diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index c201f8a1..08d3a9c8 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -555,13 +555,10 @@ def align_mumap(
                 regdct = nimpa.coreg_spm(fpet, fute,
                                          outpath=os.path.join(outpath, 'PET', 'positioning'))
             elif reg_tool == 'niftyreg':
-                if not os.path.exists(Cnt['REGPATH']):
-                    raise ValueError('e> no valid NiftyReg executable')
                 regdct = nimpa.affine_niftyreg(
                     fpet,
                     fute,
                     outpath=os.path.join(outpath, 'PET', 'positioning'),
-                    executable=Cnt['REGPATH'],
                     omp=multiprocessing.cpu_count() / 2,                 # pcomment=fcomment,
                     rigOnly=True,
                     affDirect=False,
@@ -600,13 +597,10 @@ def align_mumap(
                 regdct = nimpa.coreg_spm(fpet, ft1w,
                                          outpath=os.path.join(outpath, 'PET', 'positioning'))
             elif reg_tool == 'niftyreg':
-                if not os.path.exists(Cnt['REGPATH']):
-                    raise ValueError('e> no valid NiftyReg executable')
                 regdct = nimpa.affine_niftyreg(
                     fpet,
                     ft1w,
                     outpath=os.path.join(outpath, 'PET', 'positioning'),
-                    executable=Cnt['REGPATH'],
                     omp=multiprocessing.cpu_count() / 2,
                     rigOnly=True,
                     affDirect=False,
@@ -682,8 +676,7 @@ def align_mumap(
     elif reg_tool == 'dipy':
         nimpa.resample_dipy(fpet, fflo, faff=faff_mrpet, fimout=freg)
     else:
-        nimpa.resample_niftyreg(fpet, fflo, faff_mrpet, fimout=freg, executable=Cnt['RESPATH'],
-                                verbose=verbose)
+        nimpa.resample_niftyreg(fpet, fflo, faff_mrpet, fimout=freg, verbose=verbose)
 
     # -get the NIfTI of registered image
     nim = nib.load(freg)
diff --git a/niftypet/nipet/img/pipe.py b/niftypet/nipet/img/pipe.py
index 351d4206..46f164b6 100644
--- a/niftypet/nipet/img/pipe.py
+++ b/niftypet/nipet/img/pipe.py
@@ -2,7 +2,6 @@
 import logging
 import os
 from numbers import Integral
-from subprocess import call
 from textwrap import dedent
 
 import numpy as np
@@ -330,16 +329,8 @@ def mmrchain(
             nimpa.create_dir(fmureg)
             # the converted nii image resample to the reference size
             fmu = os.path.join(fmureg, 'mumap_dyn_frm' + str(ifrm) + fcomment + '.nii.gz')
-            # command for resampling
-            if os.path.isfile(Cnt['RESPATH']):
-                cmd = [
-                    Cnt['RESPATH'], '-ref', fmuref, '-flo', muod['fim'], '-trans', faff_frms[ifrm],
-                    '-res', fmu, '-pad', '0']
-                if log.getEffectiveLevel() > log.INFO:
-                    cmd.append('-voff')
-                call(cmd)
-            else:
-                raise IOError('Incorrect path to NiftyReg (resampling) executable.')
+            # TODO: add `-pad 0`, drop `-inter 1`?
+            nimpa.resample_niftyreg(fmuref, muod['fim'], faff_frms[ifrm], fimout=fmu)
             # get the new mu-map from the just resampled file
             muodct = nimpa.getnii(fmu, output='all')
             muo = muodct['im']
diff --git a/setup.cfg b/setup.cfg
index 79333782..919a1371 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -50,7 +50,7 @@ install_requires=
     cuvec>=2.10.0
     miutil>=0.6.0
     nibabel>=2.4.0
-    nimpa>=2.0.0
+    nimpa>=2.4.0
     ninst>=0.7.0
     numpy>=1.14
     pydicom>=1.0.2

From 237c4d12c56f8385421e132f71b9dab5d725449f Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Sun, 21 Aug 2022 04:54:53 +0100
Subject: [PATCH 27/32] fix external PRs

---
 .github/workflows/test.yml | 4 ++--
 niftypet/__init__.py       | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index 7aee93e2..86bdd7c7 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -2,7 +2,7 @@ name: Test
 on: [push, pull_request]
 jobs:
   check:
-    if: github.event_name != 'pull_request' || github.repository_owner != 'NiftyPET'
+    if: github.event_name != 'pull_request' || !contains('OWNER,MEMBER,COLLABORATOR', github.event.pull_request.author_association)
     runs-on: ubuntu-latest
     name: Check
     steps:
@@ -34,7 +34,7 @@ jobs:
         EVENT: ${{ github.event_name }}
     - run: pre-commit run -a --show-diff-on-failure
   test:
-    if: github.event_name != 'pull_request' || github.repository_owner != 'NiftyPET'
+    if: github.event_name != 'pull_request' || !contains('OWNER,MEMBER,COLLABORATOR', github.event.pull_request.author_association)
     name: Test py${{ matrix.python }}
     runs-on: [self-hosted, python, cuda, matlab]
     strategy:
diff --git a/niftypet/__init__.py b/niftypet/__init__.py
index 8db66d3d..69e3be50 100644
--- a/niftypet/__init__.py
+++ b/niftypet/__init__.py
@@ -1 +1 @@
-__path__ = __import__("pkgutil").extend_path(__path__, __name__)
+__path__ = __import__('pkgutil').extend_path(__path__, __name__)

From 85ed18725760d129769e2d5f704ad085dcd457e9 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Thu, 15 Sep 2022 13:10:35 +0100
Subject: [PATCH 28/32] build/test: misc updates, bump ninst, add py3.10

---
 .github/workflows/comment-bot.yml | 22 ++++++++++------------
 .github/workflows/test.yml        | 16 ++++++++--------
 .pre-commit-config.yaml           |  2 +-
 pyproject.toml                    |  2 +-
 setup.cfg                         |  5 +++--
 5 files changed, 23 insertions(+), 24 deletions(-)

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/.github/workflows/test.yml b/.github/workflows/test.yml
index 86bdd7c7..9e1d7ce7 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -6,11 +6,11 @@ jobs:
     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') }}
@@ -35,13 +35,13 @@ jobs:
     - run: pre-commit run -a --show-diff-on-failure
   test:
     if: github.event_name != 'pull_request' || !contains('OWNER,MEMBER,COLLABORATOR', github.event.pull_request.author_association)
-    name: Test py${{ matrix.python }}
     runs-on: [self-hosted, python, cuda, matlab]
     strategy:
       matrix:
-        python: [3.7, 3.9]
+        python: [3.7, '3.10']
+    name: Test py${{ matrix.python }}
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@v3
       with:
         fetch-depth: 0
     - name: Run setup-python
@@ -57,10 +57,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:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 929980db..b157f736 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.2.0
+  rev: v4.3.0
   hooks:
   - id: check-added-large-files
   - id: check-case-conflict
diff --git a/pyproject.toml b/pyproject.toml
index 93efe209..554a3f78 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [build-system]
 requires = ["setuptools>=42", "wheel", "setuptools_scm[toml]>=3.4",
-            "cuvec>=2.8.0", "ninst>=0.10.0", "numpy>=1.14", "miutil[cuda]>=0.4.0",
+            "cuvec>=2.8.0", "ninst>=0.12.0", "numpy>=1.14", "miutil[cuda]>=0.4.0",
             "scikit-build>=0.11.0", "cmake>=3.18", "ninja"]
 
 [tool.setuptools_scm]
diff --git a/setup.cfg b/setup.cfg
index 919a1371..8b38d2e8 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -31,6 +31,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.
 [options]
@@ -41,7 +42,7 @@ setup_requires=
     setuptools_scm[toml]
     cuvec>=2.8.0
     miutil[cuda]>=0.4.0
-    ninst>=0.10.0
+    ninst>=0.12.0
     numpy>=1.14
     scikit-build>=0.11.0
     cmake>=3.18
@@ -51,7 +52,7 @@ install_requires=
     miutil>=0.6.0
     nibabel>=2.4.0
     nimpa>=2.4.0
-    ninst>=0.7.0
+    ninst>=0.12.0
     numpy>=1.14
     pydicom>=1.0.2
     setuptools

From 2e1ad588883dc10795a57e1c282db825d1d8a980 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Thu, 15 Sep 2022 13:15:23 +0100
Subject: [PATCH 29/32] lint

---
 niftypet/nipet/img/mmrimg.py | 36 +++++++++++++++++++++++++++++++++---
 1 file changed, 33 insertions(+), 3 deletions(-)

diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index 08d3a9c8..e3235d9b 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -559,7 +559,7 @@ def align_mumap(
                     fpet,
                     fute,
                     outpath=os.path.join(outpath, 'PET', 'positioning'),
-                    omp=multiprocessing.cpu_count() / 2,                 # pcomment=fcomment,
+                    omp=multiprocessing.cpu_count() // 2,                # pcomment=fcomment,
                     rigOnly=True,
                     affDirect=False,
                     maxit=5,
@@ -601,7 +601,7 @@ def align_mumap(
                     fpet,
                     ft1w,
                     outpath=os.path.join(outpath, 'PET', 'positioning'),
-                    omp=multiprocessing.cpu_count() / 2,
+                    omp=multiprocessing.cpu_count() // 2,
                     rigOnly=True,
                     affDirect=False,
                     maxit=5,
@@ -1108,6 +1108,36 @@ def rmumaps(datain, Cnt, t0=0, t1=0, use_stored=False):
         faff = os.path.join(os.path.dirname(ft1w), fcomment + 'mr2pet_affine' + '.txt')
         # time.strftime('%d%b%y_%H.%M',time.gmtime())
         # > call the registration routine
+
+        # TODO: add `-pad 0`, drop `-inter 1`?
+        nimpa.affine_niftyreg(
+            recute.fpet,
+            ft1w,
+            rigOnly=True,
+            speed=True,
+            outpath=os.path.dirname(ft1w),
+            fname_aff=fcomment + 'mr2pet_affine' + '.txt',
+            pi=pi,
+            pv=pv,
+            smof=smof,
+            smor=smor,
+            maxit=10,
+            omp=multiprocessing.cpu_count() // 2,
+        )
+        # pickname='ref',
+        # fcomment='',
+        # executable=None,
+        # omp=1,
+        # affDirect=False,
+        # maxit=5,
+        # rmsk=True,
+        # fmsk=True,
+        # rfwhm=15.,
+        # rthrsh=0.05,
+        # ffwhm=15.,
+        # fthrsh=0.05,
+        # verbose=True,
+
         if os.path.isfile(Cnt['REGPATH']):
             cmd = [
                 Cnt['REGPATH'], '-ref', recute.fpet, '-flo', ft1w, '-rigOnly', '-speeeeed', '-aff',
@@ -1229,7 +1259,7 @@ def rmumaps(datain, Cnt, t0=0, t1=0, use_stored=False):
 #                 ft1w,
 #                 outpath=os.path.join(outpath, 'PET', 'positioning'), # pcomment=fcomment,
 #                 executable=Cnt['REGPATH'],
-#                 omp=multiprocessing.cpu_count() / 2,
+#                 omp=multiprocessing.cpu_count() // 2,
 #                 rigOnly=True,
 #                 affDirect=False,
 #                 maxit=5,

From 9365f59001dee78d33d09a0d22154b37051de327 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Thu, 15 Sep 2022 13:25:46 +0100
Subject: [PATCH 30/32] scipy.ndimage.filters => scipy.ndimage

---
 examples/demo.ipynb          |  2 +-
 niftypet/nipet/lm/mmrhist.py |  2 +-
 niftypet/nipet/prj/mmrrec.py |  4 ++--
 niftypet/nipet/sct/mmrsct.py | 14 +++++++-------
 4 files changed, 11 insertions(+), 11 deletions(-)

diff --git a/examples/demo.ipynb b/examples/demo.ipynb
index a03def75..438d01c6 100644
--- a/examples/demo.ipynb
+++ b/examples/demo.ipynb
@@ -51,7 +51,7 @@
     "import cuvec as cu\n",
     "import matplotlib.pyplot as plt\n",
     "import numpy as np\n",
-    "from scipy.ndimage.filters import gaussian_filter\n",
+    "from scipy.ndimage import gaussian_filter\n",
     "from tqdm.auto import trange\n",
     "\n",
     "from niftypet import nipet, nimpa\n",
diff --git a/niftypet/nipet/lm/mmrhist.py b/niftypet/nipet/lm/mmrhist.py
index 20a99f3d..b347c058 100644
--- a/niftypet/nipet/lm/mmrhist.py
+++ b/niftypet/nipet/lm/mmrhist.py
@@ -132,7 +132,7 @@ def hist(
     pvs_sgtl = np.right_shift(hstout['pvs'], 8).astype(np.float32)
     pvs_crnl = np.bitwise_and(hstout['pvs'], 255).astype(np.float32)
 
-    cmass = Cnt['SO_VXZ'] * ndi.filters.gaussian_filter(hstout['mss'], cmass_sig, mode='mirror')
+    cmass = Cnt['SO_VXZ'] * ndi.gaussian_filter(hstout['mss'], cmass_sig, mode='mirror')
     log.debug(
         'centre of mass of axial radiodistribution (filtered with Gaussian of SD ={}):  COMPLETED.'
         .format(cmass_sig))
diff --git a/niftypet/nipet/prj/mmrrec.py b/niftypet/nipet/prj/mmrrec.py
index 62cd7f99..a92b16eb 100644
--- a/niftypet/nipet/prj/mmrrec.py
+++ b/niftypet/nipet/prj/mmrrec.py
@@ -431,8 +431,8 @@ def osemone(datain, mumaps, hst, scanner_params, recmod=3, itr=4, fwhm=0., psf=N
     im_smo = None
     fsmo = None
     if fwhm > 0:
-        im_smo = ndi.filters.gaussian_filter(im, fwhm2sig(fwhm, voxsize=Cnt['SZ_VOXY'] * 10),
-                                             mode='mirror')
+        im_smo = ndi.gaussian_filter(im, fwhm2sig(fwhm, voxsize=Cnt['SZ_VOXY'] * 10),
+                                     mode='mirror')
 
         if store_img:
             fsmo = fpet.split('.nii.gz')[0] + '_smo-' + str(fwhm).replace('.', '-') + 'mm.nii.gz'
diff --git a/niftypet/nipet/sct/mmrsct.py b/niftypet/nipet/sct/mmrsct.py
index de830b48..8f5dfebf 100644
--- a/niftypet/nipet/sct/mmrsct.py
+++ b/niftypet/nipet/sct/mmrsct.py
@@ -445,8 +445,8 @@ def vsm(
 
     # -smooth for defining the sino scatter only regions
     if fwhm_input > 0.:
-        mu_sctonly = ndi.filters.gaussian_filter(mmrimg.convert2dev(muo, Cnt),
-                                                 fwhm2sig(fwhm_input, Cnt), mode='mirror')
+        mu_sctonly = ndi.gaussian_filter(mmrimg.convert2dev(muo, Cnt), fwhm2sig(fwhm_input, Cnt),
+                                         mode='mirror')
     else:
         mu_sctonly = muo
 
@@ -466,8 +466,8 @@ def vsm(
 
     # > smooth before scaling/down-sampling the mu-map and emission images
     if fwhm_input > 0.:
-        muim = ndi.filters.gaussian_filter(muo + muh, fwhm2sig(fwhm_input, Cnt), mode='mirror')
-        emim = ndi.filters.gaussian_filter(em, fwhm2sig(fwhm_input, Cnt), mode='mirror')
+        muim = ndi.gaussian_filter(muo + muh, fwhm2sig(fwhm_input, Cnt), mode='mirror')
+        emim = ndi.gaussian_filter(em, fwhm2sig(fwhm_input, Cnt), mode='mirror')
     else:
         muim = muo + muh
         emim = em
@@ -478,7 +478,7 @@ def vsm(
     # -smooth the mu-map for mask creation.
     # the mask contains voxels for which attenuation ray LUT is found.
     if fwhm_input > 0.:
-        smomu = ndi.filters.gaussian_filter(muim, fwhm2sig(fwhm_input, Cnt), mode='mirror')
+        smomu = ndi.gaussian_filter(muim, fwhm2sig(fwhm_input, Cnt), mode='mirror')
         mumsk = np.int8(smomu > 0.003)
     else:
         mumsk = np.int8(muim > 0.001)
@@ -608,9 +608,9 @@ def vsm(
         eim = eim[:, ::-1, ::-1]
         eim = np.transpose(eim, (2, 1, 0))
 
-        em_sctonly = ndi.filters.gaussian_filter(eim, fwhm2sig(.6, Cnt), mode='mirror')
+        em_sctonly = ndi.gaussian_filter(eim, fwhm2sig(.6, Cnt), mode='mirror')
         msk = np.float32(em_sctonly > 0.07 * np.max(em_sctonly))
-        msk = ndi.filters.gaussian_filter(msk, fwhm2sig(.6, Cnt), mode='mirror')
+        msk = ndi.gaussian_filter(msk, fwhm2sig(.6, Cnt), mode='mirror')
         msk = np.float32(msk > 0.01)
         msksn = mmrprj.frwd_prj(msk, txLUT, axLUT, Cnt)
 

From 766a1e0f4f79a1f6afb62eac28d8d868340f4efc Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Thu, 27 Oct 2022 23:35:07 +0100
Subject: [PATCH 31/32] fixup nimpa.affine_niftyreg

---
 niftypet/nipet/img/mmrimg.py | 58 +++++++++---------------------------
 1 file changed, 14 insertions(+), 44 deletions(-)

diff --git a/niftypet/nipet/img/mmrimg.py b/niftypet/nipet/img/mmrimg.py
index e3235d9b..ed3e0c4a 100644
--- a/niftypet/nipet/img/mmrimg.py
+++ b/niftypet/nipet/img/mmrimg.py
@@ -1,5 +1,4 @@
 """Image functions for PET data reconstruction and processing."""
-
 import glob
 import logging
 import math
@@ -7,7 +6,6 @@
 import os
 import re
 import shutil
-from subprocess import run
 
 import nibabel as nib
 import numpy as np
@@ -1102,54 +1100,26 @@ def rmumaps(datain, Cnt, t0=0, t1=0, use_stored=False):
         else:
             raise IOError('Disaster: no T1w image!')
 
-        # putput for the T1w in register with PET
-        ft1out = os.path.join(os.path.dirname(ft1w), 'T1w_r' + '.nii.gz')
-        # pext file fo rthe affine transform T1w->PET
-        faff = os.path.join(os.path.dirname(ft1w), fcomment + 'mr2pet_affine' + '.txt')
+        # > putput for the T1w in register with PET
+        # ft1out = os.path.join(os.path.dirname(ft1w), 'T1w_r' + '.nii.gz')
+        # > pext file fo rthe affine transform T1w->PET
+        # faff = os.path.join(os.path.dirname(ft1w), fcomment + 'mr2pet_affine' + '.txt')
         # time.strftime('%d%b%y_%H.%M',time.gmtime())
         # > call the registration routine
-
+        # cmd = [
+        #     Cnt['REGPATH'], '-ref', recute.fpet, '-flo', ft1w, '-rigOnly', '-speeeeed', '-aff',
+        #     faff, '-res', ft1out]
         # TODO: add `-pad 0`, drop `-inter 1`?
-        nimpa.affine_niftyreg(
-            recute.fpet,
-            ft1w,
-            rigOnly=True,
-            speed=True,
-            outpath=os.path.dirname(ft1w),
-            fname_aff=fcomment + 'mr2pet_affine' + '.txt',
-            pi=pi,
-            pv=pv,
-            smof=smof,
-            smor=smor,
-            maxit=10,
-            omp=multiprocessing.cpu_count() // 2,
-        )
-        # pickname='ref',
-        # fcomment='',
-        # executable=None,
-        # omp=1,
-        # affDirect=False,
+        # pi=50, pv=50, smof=0, smor=0,
         # maxit=5,
-        # rmsk=True,
-        # fmsk=True,
-        # rfwhm=15.,
-        # rthrsh=0.05,
-        # ffwhm=15.,
-        # fthrsh=0.05,
-        # verbose=True,
-
-        if os.path.isfile(Cnt['REGPATH']):
-            cmd = [
-                Cnt['REGPATH'], '-ref', recute.fpet, '-flo', ft1w, '-rigOnly', '-speeeeed', '-aff',
-                faff, '-res', ft1out]
-            if log.getEffectiveLevel() > logging.INFO:
-                cmd.append('-voff')
-            run(cmd)
-        else:
-            raise IOError('Path to registration executable is incorrect!')
+        regdct = nimpa.affine_niftyreg(recute.fpet, ft1w, rigOnly=True, speed=True,
+                                       outpath=os.path.dirname(ft1w),
+                                       fname_aff=fcomment + 'mr2pet_affine' + '.txt',
+                                       omp=multiprocessing.cpu_count() // 2,
+                                       verbose=log.getEffectiveLevel() <= logging.INFO)
 
         # pet the pCT mu-map with the above faff
-        pmudic = pct_mumap(datain, txLUT_, axLUT_, Cnt, faff=faff, fpet=recute.fpet,
+        pmudic = pct_mumap(datain, txLUT_, axLUT_, Cnt, faff=regdct['faff'], fpet=recute.fpet,
                            fcomment=fcomment)
         mup = pmudic['im']
 

From eb8e98aef02115139fe551f3acf2592fe5eaee61 Mon Sep 17 00:00:00 2001
From: Casper da Costa-Luis <imaging@cdcl.ml>
Date: Thu, 27 Oct 2022 23:35:33 +0100
Subject: [PATCH 32/32] minor readme installation update

---
 README.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.rst b/README.rst
index e51061e6..f8b58bb1 100644
--- a/README.rst
+++ b/README.rst
@@ -37,7 +37,7 @@ It's also recommended (but not required) to use `conda`.
     export HMUDIR=$HOME/mmr_hardwareumaps
     # cross-platform install
     conda install -c conda-forge python=3 \
-      ipykernel numpy scipy scikit-image matplotlib ipywidgets
+      ipykernel numpy scipy scikit-image matplotlib ipywidgets dipy nibabel pydicom
     pip install "nipet>=2"
 
 External CMake Projects