diff --git a/Dockerfile b/Dockerfile
new file mode 100644
index 0000000..2e0aa4b
--- /dev/null
+++ b/Dockerfile
@@ -0,0 +1,9 @@
+FROM python:3.8-slim
+
+COPY . /gradunwarp
+
+RUN cd /gradunwarp \
+  && pip install -r requirements.txt \
+  && pip install .
+
+ENTRYPOINT ["/usr/local/bin/gradient_unwarp.py"]
diff --git a/gradunwarp/core/coeffs.py b/gradunwarp/core/coeffs.py
index efa9a77..9e1f8ca 100644
--- a/gradunwarp/core/coeffs.py
+++ b/gradunwarp/core/coeffs.py
@@ -13,46 +13,52 @@
 from .globals import siemens_cas, ge_cas
 
 
-log = logging.getLogger('gradunwarp')
+log = logging.getLogger("gradunwarp")
 
 
-Coeffs = namedtuple('Coeffs', 'alpha_x, alpha_y, alpha_z, \
-                        beta_x, beta_y, beta_z, R0_m')
+Coeffs = namedtuple(
+    "Coeffs",
+    "alpha_x, alpha_y, alpha_z, \
+                        beta_x, beta_y, beta_z, R0_m",
+)
 
 try:
     advance_iterator = next
 except NameError:
+
     def advance_iterator(it):
         return it.next()
+
+
 next = advance_iterator
 
 
 def get_coefficients(vendor, cfile):
-    ''' depending on the vendor and the coefficient file,
+    """depending on the vendor and the coefficient file,
     return the spherical harmonics coefficients as a named tuple.
-    '''
-    log.info('Parsing ' + cfile + ' for harmonics coeffs')
-    if vendor == 'siemens' and cfile.endswith('.coef'):
+    """
+    log.info("Parsing " + cfile + " for harmonics coeffs")
+    if vendor == "siemens" and cfile.endswith(".coef"):
         return get_siemens_coef(cfile)
-    if vendor == 'siemens' and cfile.endswith('.grad'):
+    if vendor == "siemens" and cfile.endswith(".grad"):
         return get_siemens_grad(cfile)
 
 
 def coef_file_parse(cfile, txt_var_map):
-    ''' a separate function because GE and Siemens .coef files
+    """a separate function because GE and Siemens .coef files
     have similar structure
 
     modifies txt_var_map in place
-    '''
+    """
     # parse .coef file. Strip unneeded characters. a valid line in that file is
     # broken into validline_list
-    coef_re = re.compile('^[^\#]')  # regex for first character not a '#'
-    coef_file = open(cfile, 'r')
+    coef_re = re.compile("^[^\#]")  # regex for first character not a '#'
+    coef_file = open(cfile, "r")
     for line in coef_file.readlines():
         if coef_re.match(line):
-            validline_list = line.lstrip(' \t').rstrip(';\n').split()
+            validline_list = line.lstrip(" \t").rstrip(";\n").split()
             if validline_list:
-                log.info('Parsed : %s' % validline_list)
+                log.info("Parsed : %s" % validline_list)
                 l = validline_list
                 x = int(l[1])
                 y = int(l[2])
@@ -60,23 +66,25 @@ def coef_file_parse(cfile, txt_var_map):
 
 
 def get_siemens_coef(cfile):
-    ''' Parse the Siemens .coef file.
+    """Parse the Siemens .coef file.
     Note that R0_m is not explicitly contained in the file
-    '''
-    R0m_map = {'sonata': 0.25,
-               'avanto': 0.25,
-               'quantum': 0.25,
-               'allegra': 0.14,
-               'as39s': 0.25,
-               'as39st': 0.25,
-               'as39t': 0.25}
+    """
+    R0m_map = {
+        "sonata": 0.25,
+        "avanto": 0.25,
+        "quantum": 0.25,
+        "allegra": 0.14,
+        "as39s": 0.25,
+        "as39st": 0.25,
+        "as39t": 0.25,
+    }
     for rad in R0m_map.keys():
         if cfile.startswith(rad):
             R0_m = R0m_map[rad]
 
     coef_array_sz = siemens_cas
     # allegra is slightly different
-    if cfile.startswith('allegra'):
+    if cfile.startswith("allegra"):
         coef_array_sz = 15
     ax = np.zeros((coef_array_sz, coef_array_sz))
     ay = np.zeros((coef_array_sz, coef_array_sz))
@@ -84,12 +92,14 @@ def get_siemens_coef(cfile):
     bx = np.zeros((coef_array_sz, coef_array_sz))
     by = np.zeros((coef_array_sz, coef_array_sz))
     bz = np.zeros((coef_array_sz, coef_array_sz))
-    txt_var_map = {'Alpha_x': ax,
-                   'Alpha_y': ay,
-                   'Alpha_z': az,
-                   'Beta_x': bx,
-                   'Beta_y': by,
-                   'Beta_z': bz}
+    txt_var_map = {
+        "Alpha_x": ax,
+        "Alpha_y": ay,
+        "Alpha_z": az,
+        "Beta_x": bx,
+        "Beta_y": by,
+        "Beta_z": bz,
+    }
 
     coef_file_parse(cfile, txt_var_map)
 
@@ -97,35 +107,37 @@ def get_siemens_coef(cfile):
 
 
 def get_ge_coef(cfile):
-    ''' Parse the GE .coef file.
-    '''
+    """Parse the GE .coef file."""
     ax = np.zeros((ge_cas, ge_cas))
     ay = np.zeros((ge_cas, ge_cas))
     az = np.zeros((ge_cas, ge_cas))
     bx = np.zeros((ge_cas, ge_cas))
     by = np.zeros((ge_cas, ge_cas))
     bz = np.zeros((ge_cas, ge_cas))
-    txt_var_map = {'Alpha_x': ax,
-                   'Alpha_y': ay,
-                   'Alpha_z': az,
-                   'Beta_x': bx,
-                   'Beta_y': by,
-                   'Beta_z': bz}
+    txt_var_map = {
+        "Alpha_x": ax,
+        "Alpha_y": ay,
+        "Alpha_z": az,
+        "Beta_x": bx,
+        "Beta_y": by,
+        "Beta_z": bz,
+    }
 
     coef_file_parse(cfile, txt_var_map)
 
     return Coeffs(ax, ay, az, bx, by, bz, R0_m)
 
+
 def grad_file_parse(gfile, txt_var_map):
-    ''' a separate function because GE and Siemens .coef files
+    """a separate function because GE and Siemens .coef files
     have similar structure
 
     modifies txt_var_map in place
-    '''
-    gf = open(gfile, 'r')
+    """
+    gf = open(gfile, "r")
     line = next(gf)
     # skip the comments
-    while not line.startswith('#*] END:'):
+    while not line.startswith("#*] END:"):
         line = next(gf)
 
     # get R0
@@ -147,30 +159,30 @@ def grad_file_parse(gfile, txt_var_map):
     ymax = 0
 
     while 1:
-        lindex =  line.find('(')
-        rindex =  line.find(')')
+        lindex = line.find("(")
+        rindex = line.find(")")
         if lindex == -1 and rindex == -1:
             break
-        arrindex = line[lindex+1:rindex]
-        xs, ys = arrindex.split(',')
-        x = int(xs) 
+        arrindex = line[lindex + 1 : rindex]
+        xs, ys = arrindex.split(",")
+        x = int(xs)
         y = int(ys)
         if x > xmax:
             xmax = x
         if y > ymax:
             ymax = y
-        if line.find('A') != -1 and line.find('x') != -1:
-            txt_var_map['Alpha_x'][x,y] = float(line.split()[-2])
-        if line.find('A') != -1 and line.find('y') != -1:
-            txt_var_map['Alpha_y'][x,y] = float(line.split()[-2])
-        if line.find('A') != -1 and line.find('z') != -1:
-            txt_var_map['Alpha_z'][x,y] = float(line.split()[-2])
-        if line.find('B') != -1 and line.find('x') != -1:
-            txt_var_map['Beta_x'][x,y] = float(line.split()[-2])
-        if line.find('B') != -1 and line.find('y') != -1:
-            txt_var_map['Beta_y'][x,y] = float(line.split()[-2])
-        if line.find('B') != -1 and line.find('z') != -1:
-            txt_var_map['Beta_z'][x,y] = float(line.split()[-2])
+        if line.find("A") != -1 and line.find("x") != -1:
+            txt_var_map["Alpha_x"][x, y] = float(line.split()[-2])
+        if line.find("A") != -1 and line.find("y") != -1:
+            txt_var_map["Alpha_y"][x, y] = float(line.split()[-2])
+        if line.find("A") != -1 and line.find("z") != -1:
+            txt_var_map["Alpha_z"][x, y] = float(line.split()[-2])
+        if line.find("B") != -1 and line.find("x") != -1:
+            txt_var_map["Beta_x"][x, y] = float(line.split()[-2])
+        if line.find("B") != -1 and line.find("y") != -1:
+            txt_var_map["Beta_y"][x, y] = float(line.split()[-2])
+        if line.find("B") != -1 and line.find("z") != -1:
+            txt_var_map["Beta_z"][x, y] = float(line.split()[-2])
         try:
             line = next(gf)
         except StopIteration:
@@ -179,12 +191,12 @@ def grad_file_parse(gfile, txt_var_map):
     # just return R0_m but also txt_var_map is returned
     return R0_m, (xmax, ymax)
 
+
 def get_siemens_grad(gfile):
-    ''' Parse the siemens .grad file
-    '''
+    """Parse the siemens .grad file"""
     coef_array_sz = siemens_cas
     # allegra is slightly different
-    if gfile.startswith('coef_AC44'):
+    if gfile.startswith("coef_AC44"):
         coef_array_sz = 15
     ax = np.zeros((coef_array_sz, coef_array_sz))
     ay = np.zeros((coef_array_sz, coef_array_sz))
@@ -192,23 +204,24 @@ def get_siemens_grad(gfile):
     bx = np.zeros((coef_array_sz, coef_array_sz))
     by = np.zeros((coef_array_sz, coef_array_sz))
     bz = np.zeros((coef_array_sz, coef_array_sz))
-    txt_var_map = {'Alpha_x': ax,
-                   'Alpha_y': ay,
-                   'Alpha_z': az,
-                   'Beta_x': bx,
-                   'Beta_y': by,
-                   'Beta_z': bz}
+    txt_var_map = {
+        "Alpha_x": ax,
+        "Alpha_y": ay,
+        "Alpha_z": az,
+        "Beta_x": bx,
+        "Beta_y": by,
+        "Beta_z": bz,
+    }
 
     R0_m, max_ind = grad_file_parse(gfile, txt_var_map)
     ind = max(max_ind)
 
     # pruned alphas and betas
-    ax = ax[:ind+1, :ind+1]
-    ay = ay[:ind+1, :ind+1]
-    az = az[:ind+1, :ind+1]
-    bx = bx[:ind+1, :ind+1]
-    by = by[:ind+1, :ind+1]
-    bz = bz[:ind+1, :ind+1]
+    ax = ax[: ind + 1, : ind + 1]
+    ay = ay[: ind + 1, : ind + 1]
+    az = az[: ind + 1, : ind + 1]
+    bx = bx[: ind + 1, : ind + 1]
+    by = by[: ind + 1, : ind + 1]
+    bz = bz[: ind + 1, : ind + 1]
 
     return Coeffs(ax, ay, az, bx, by, bz, R0_m)
-
diff --git a/gradunwarp/core/globals.py b/gradunwarp/core/globals.py
index a07a4e0..68d4061 100644
--- a/gradunwarp/core/globals.py
+++ b/gradunwarp/core/globals.py
@@ -6,20 +6,20 @@
 ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
 import logging
 
-VERSION = '1.2.1+HCP'
+VERSION = "1.2.1+HCP"
 
-usage = '''
+usage = """
 gradient_unwarp infile outfile manufacturer -g <coefficient file> [optional arguments]
-'''
+"""
 
 
 # SIEMENS stuff
 siemens_cas = 100  # coefficient array size
-siemens_fovmin = -.30  # fov min in meters
-siemens_fovmax = .30  # fov max in meters
-siemens_numpoints = 60 # number of grid points in each direction
+siemens_fovmin = -0.30  # fov min in meters
+siemens_fovmax = 0.30  # fov max in meters
+siemens_numpoints = 60  # number of grid points in each direction
 # max jacobian determinant for siemens
-siemens_max_det = 10.
+siemens_max_det = 10.0
 
 
 # GE stuff
@@ -30,7 +30,7 @@
 
 
 def get_logger():
-    log = logging.getLogger('gradunwarp')
+    log = logging.getLogger("gradunwarp")
     log.setLevel(logging.INFO)
     ch = logging.StreamHandler()
     ch.setLevel(logging.INFO)
diff --git a/gradunwarp/core/gradient_unwarp.py b/gradunwarp/core/gradient_unwarp.py
index 49aee18..2f225fd 100755
--- a/gradunwarp/core/gradient_unwarp.py
+++ b/gradunwarp/core/gradient_unwarp.py
@@ -8,110 +8,141 @@
 import argparse as arg
 import os
 import logging
-from gradunwarp.core import (globals, coeffs, utils)
+from gradunwarp.core import globals, coeffs, utils
 from gradunwarp.core.unwarp_resample import Unwarper
+import nibabel as nb
 
 log = globals.get_logger()
 
 
 def argument_parse_gradunwarp():
-    '''Arguments parser from the command line
-    '''
+    """Arguments parser from the command line"""
     # initiate
     try:
         p = arg.ArgumentParser(usage=globals.usage)
-        p.add_argument('--version', '-v', action='version', version=globals.VERSION)
+        p.add_argument("--version", "-v", action="version", version=globals.VERSION)
     except:
-        #maintain compatibility with pre py2.7 argparse (deprecated in py2.7 but fails in py3)
+        # maintain compatibility with pre py2.7 argparse (deprecated in py2.7 but fails in py3)
         p = arg.ArgumentParser(version=globals.VERSION, usage=globals.usage)
 
     # required arguments
-    p.add_argument('infile', action='store',
-                  help='The input warped file (nifti or mgh)')
-    p.add_argument('outfile', action='store',
-                  help='The output unwarped file (extension should be .nii/.nii.gz/.mgh/.mgz)')
-    p.add_argument('vendor', action='store', choices=['siemens', 'ge'], 
-                  help='vendor (either "ge" or "siemens" for now)')
+    p.add_argument(
+        "infile", action="store", help="The input warped file (nifti or mgh)"
+    )
+    p.add_argument(
+        "outfile",
+        action="store",
+        help="The output unwarped file (extension should be .nii/.nii.gz/.mgh/.mgz)",
+    )
+    p.add_argument(
+        "vendor",
+        action="store",
+        choices=["siemens", "ge"],
+        help='vendor (either "ge" or "siemens" for now)',
+    )
 
     coef_grp = p.add_mutually_exclusive_group(required=True)
-    coef_grp.add_argument('-g', '--gradfile', dest='gradfile',
-                         help='The .grad coefficient file')
-    coef_grp.add_argument('-c', '--coeffile', dest='coeffile',
-                         help='The .coef coefficient file')
+    coef_grp.add_argument(
+        "-g", "--gradfile", dest="gradfile", help="The .grad coefficient file"
+    )
+    coef_grp.add_argument(
+        "-c", "--coeffile", dest="coeffile", help="The .coef coefficient file"
+    )
 
     # optional arguments
-    p.add_argument('-w', '--warp', action='store_true', default=False,
-                  help='warp a volume (as opposed to unwarping)')
-    p.add_argument('-n', '--nojacobian', dest='nojac', action='store_true',
-                  default=False, help='Do not perform Jacobian intensity correction')
-    p.add_argument('--fovmin', dest='fovmin',
-                  help='the minimum extent of harmonics evaluation grid in meters')
-    p.add_argument('--fovmax', dest='fovmax',
-                  help='the maximum extent of harmonics evaluation grid in meters')
-    p.add_argument('--numpoints', dest='numpoints',
-                   help='number of grid points in each direction')
-    p.add_argument('--interp_order', dest='order',
-                   help='the order of interpolation(1..4) where 1 is linear - default')
-
-    p.add_argument('--verbose', action='store_true', default=False)
+    p.add_argument(
+        "-w",
+        "--warp",
+        action="store_true",
+        default=False,
+        help="warp a volume (as opposed to unwarping)",
+    )
+    p.add_argument(
+        "-n",
+        "--nojacobian",
+        dest="nojac",
+        action="store_true",
+        default=False,
+        help="Do not perform Jacobian intensity correction",
+    )
+    p.add_argument(
+        "--fovmin",
+        dest="fovmin",
+        help="the minimum extent of harmonics evaluation grid in meters",
+    )
+    p.add_argument(
+        "--fovmax",
+        dest="fovmax",
+        help="the maximum extent of harmonics evaluation grid in meters",
+    )
+    p.add_argument(
+        "--numpoints", dest="numpoints", help="number of grid points in each direction"
+    )
+    p.add_argument(
+        "--interp_order",
+        dest="order",
+        help="the order of interpolation(1..4) where 1 is linear - default",
+    )
+
+    p.add_argument("--verbose", action="store_true", default=False)
 
     args = p.parse_args()
 
     # do some validation
     if not os.path.exists(args.infile):
-        raise IOError(args.infile + ' not found')
+        raise IOError(args.infile + " not found")
     if args.gradfile:
         if not os.path.exists(args.gradfile):
-            raise IOError(args.gradfile + ' not found')
+            raise IOError(args.gradfile + " not found")
     if args.coeffile:
         if not os.path.exists(args.coeffile):
-            raise IOError(args.coeffile + ' not found')
+            raise IOError(args.coeffile + " not found")
 
     return args
 
 
 class GradientUnwarpRunner(object):
-    ''' Takes the option datastructure after parsing the commandline.
+    """Takes the option datastructure after parsing the commandline.
     run() method performs the actual unwarping
     write() method performs the writing of the unwarped volume
-    '''
+    """
+
     def __init__(self, args):
-        ''' constructor takes the option datastructure which is the
+        """constructor takes the option datastructure which is the
         result of (options, args) = parser.parse_args()
-        '''
+        """
         self.args = args
         self.unwarper = None
 
         log.setLevel(logging.INFO)
-        if hasattr(self.args, 'verbose'):
+        if hasattr(self.args, "verbose"):
             log.setLevel(logging.DEBUG)
 
     def run(self):
-        ''' run the unwarp resample
-        '''
+        """run the unwarp resample"""
         # get the spherical harmonics coefficients from parsing
         # the given .coeff file xor .grad file
-        if hasattr(self.args, 'gradfile') and self.args.gradfile:
-            self.coeffs = coeffs.get_coefficients(self.args.vendor,
-                                                 self.args.gradfile)
+        if hasattr(self.args, "gradfile") and self.args.gradfile:
+            self.coeffs = coeffs.get_coefficients(self.args.vendor, self.args.gradfile)
         else:
-            self.coeffs = coeffs.get_coefficients(self.args.vendor,
-                                                 self.args.coeffile)
+            self.coeffs = coeffs.get_coefficients(self.args.vendor, self.args.coeffile)
 
-        self.vol, self.m_rcs2ras = utils.get_vol_affine(self.args.infile)
+        self.input_nii = nb.load(self.args.infile)
 
-        self.unwarper = Unwarper(self.vol, self.m_rcs2ras, self.args.vendor, self.coeffs, self.args.infile )
-        if hasattr(self.args, 'fovmin') and self.args.fovmin:
+        self.unwarper = Unwarper(
+            self.input_nii, self.args.vendor, self.coeffs, self.args.infile
+        )
+        if hasattr(self.args, "fovmin") and self.args.fovmin:
             self.unwarper.fovmin = float(self.args.fovmin)
-        if hasattr(self.args, 'fovmax') and self.args.fovmax:
+        if hasattr(self.args, "fovmax") and self.args.fovmax:
             self.unwarper.fovmax = float(self.args.fovmax)
-        if hasattr(self.args, 'numpoints') and self.args.numpoints:
+        if hasattr(self.args, "numpoints") and self.args.numpoints:
             self.unwarper.numpoints = int(self.args.numpoints)
-        if hasattr(self.args, 'warp') and self.args.warp:
+        if hasattr(self.args, "warp") and self.args.warp:
             self.unwarper.warp = True
-        if hasattr(self.args, 'nojac') and self.args.nojac:
+        if hasattr(self.args, "nojac") and self.args.nojac:
             self.unwarper.nojac = True
-        if hasattr(self.args, 'order') and self.args.order:
+        if hasattr(self.args, "order") and self.args.order:
             self.unwarper.order = int(self.args.order)
         self.unwarper.run()
 
@@ -119,7 +150,7 @@ def write(self):
         self.unwarper.write(self.args.outfile)
 
 
-if __name__ == '__main__':
+if __name__ == "__main__":
     args = argument_parse_gradunwarp()
 
     grad_unwarp = GradientUnwarpRunner(args)
diff --git a/gradunwarp/core/interp3_ext.c b/gradunwarp/core/interp3_ext.c
deleted file mode 100644
index 9990e3c..0000000
--- a/gradunwarp/core/interp3_ext.c
+++ /dev/null
@@ -1,198 +0,0 @@
-/*
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-#
-#   See COPYING file distributed along with the gradunwarp package for the
-#   copyright and license terms.
-#
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-*/
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include <math.h>
-
-# define CUBE(x)   ((x) * (x) * (x))
-# define SQR(x)    ((x) * (x))
-
-/* compile command
- gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.7 -o interp3_ext.so interp3_ext.c
- */
-static PyObject *interp3(PyObject *self, PyObject *args);
-float TriCubic (float px, float py, float pz, float *volume, int xDim, int yDim, int zDim);
-
-// what function are exported
-static PyMethodDef tricubicmethods[] = {
-	{"_interp3", interp3, METH_VARARGS},
-	{NULL, NULL}
-};
-
-// This function is essential for an extension for Numpy created in C
-#if PY_MAJOR_VERSION >= 3
-
-static struct PyModuleDef moduledef = {
-        PyModuleDef_HEAD_INIT,
-        "interp3_ext",
-        NULL,
-        -1,
-        tricubicmethods,
-        NULL,
-        NULL,
-        NULL,
-        NULL
-};
-#define INITERROR return NULL
-
-PyObject *
-PyInit_interp3_ext(void)
-
-#else
-#define INITERROR return
-
-void
-initinterp3_ext(void)
-#endif
-{
-#if PY_MAJOR_VERSION >= 3
-    PyObject *module = PyModule_Create(&moduledef);
-#else
-    PyObject *module = Py_InitModule("interp3_ext", tricubicmethods);
-#endif
-
-    if (module == NULL)
-        INITERROR;
-
-    import_array();
-
-#if PY_MAJOR_VERSION >= 3
-    return module;
-#endif
-}
-
-// the data should be FLOAT32 and should be ensured in the wrapper 
-static PyObject *interp3(PyObject *self, PyObject *args)
-{
-	PyArrayObject *volume, *result, *C, *R, *S;
-	float *pr, *pc, *ps;
-        float *pvol, *pvc;
-        int xdim, ydim, zdim;
-
-	// We expect 4 arguments of the PyArray_Type
-	if(!PyArg_ParseTuple(args, "O!O!O!O!", 
-				&PyArray_Type, &volume,
-				&PyArray_Type, &R,
-				&PyArray_Type, &C,
-				&PyArray_Type, &S)) return NULL;
-
-	if ( NULL == volume ) return NULL;
-	if ( NULL == C ) return NULL;
-	if ( NULL == R ) return NULL;
-	if ( NULL == S ) return NULL;
-
-	// result matrix is the same size as C and is float
-	result = (PyArrayObject*) PyArray_ZEROS(PyArray_NDIM(C), C->dimensions, NPY_FLOAT, 0); 
-	// This is for reference counting ( I think )
-	PyArray_FLAGS(result) |= NPY_OWNDATA; 
-
-	// massive use of iterators to progress through the data
-	PyArrayIterObject *itr_v, *itr_r, *itr_c, *itr_s;
-    itr_v = (PyArrayIterObject *) PyArray_IterNew(result);
-    itr_r = (PyArrayIterObject *) PyArray_IterNew(R);
-    itr_c = (PyArrayIterObject *) PyArray_IterNew(C);
-    itr_s = (PyArrayIterObject *) PyArray_IterNew(S);
-    pvol = (float *)PyArray_DATA(volume);
-    xdim = PyArray_DIM(volume, 0);
-    ydim = PyArray_DIM(volume, 1);
-    zdim = PyArray_DIM(volume, 2);
-    //printf("%f\n", pvol[4*20*30 + 11*30 + 15]);
-    while(PyArray_ITER_NOTDONE(itr_v)) 
-    {
-		pvc = (float *) PyArray_ITER_DATA(itr_v);
-		pr = (float *) PyArray_ITER_DATA(itr_r);
-		pc = (float *) PyArray_ITER_DATA(itr_c);
-		ps = (float *) PyArray_ITER_DATA(itr_s);
-        // The order is weird because the tricubic code below is 
-        // for Fortran ordering. Note that the xdim changes fast in
-        // the code, whereas the rightmost dim should change fast
-        // in C multidimensional arrays.
-		*pvc = TriCubic(*ps, *pc, *pr, pvol, zdim, ydim, xdim); 
-		PyArray_ITER_NEXT(itr_v);
-		PyArray_ITER_NEXT(itr_r);
-		PyArray_ITER_NEXT(itr_c);
-		PyArray_ITER_NEXT(itr_s);
-    }
-
-	return result;
-}
-		
-/*
- * TriCubic - tri-cubic interpolation at point, p.
- *   inputs:
- *     px, py, pz - the interpolation point.
- *     volume - a pointer to the float volume data, stored in x,
- *              y, then z order (x index increasing fastest).
- *     xDim, yDim, zDim - dimensions of the array of volume data.
- *   returns:
- *     the interpolated value at p.
- *   note:
- *     rudimentary range checking is done in this function.
- */
-
-float TriCubic (float px, float py, float pz, float *volume, int xDim, int yDim, int zDim)
-{
-  int             x, y, z;
-  int    i, j, k;
-  float           dx, dy, dz;
-  float *pv;
-  float           u[4], v[4], w[4];
-  float           r[4], q[4];
-  float           vox = 0;
-  int             xyDim;
-
-  xyDim = xDim * yDim;
-
-  x = (int) px, y = (int) py, z = (int) pz;
-  // necessary evil truncating at dim-2 because tricubic needs 2 more values
-  // which is criminal near edges
-  // future work includes doing trilinear for edge cases
-  // range checking is extremely important here
-  if (x < 2 || x > xDim-3 || y < 2 || y > yDim-3 || z < 2 || z > zDim-3) 
-    return (0);
-
-  dx = px - (float) x, dy = py - (float) y, dz = pz - (float) z;
-  pv = volume + (x - 1) + (y - 1) * xDim + (z - 1) * xyDim;
-
-  /* factors for Catmull-Rom interpolation */
-
-  u[0] = -0.5 * CUBE (dx) + SQR (dx) - 0.5 * dx;
-  u[1] = 1.5 * CUBE (dx) - 2.5 * SQR (dx) + 1;
-  u[2] = -1.5 * CUBE (dx) + 2 * SQR (dx) + 0.5 * dx;
-  u[3] = 0.5 * CUBE (dx) - 0.5 * SQR (dx);
-
-  v[0] = -0.5 * CUBE (dy) + SQR (dy) - 0.5 * dy;
-  v[1] = 1.5 * CUBE (dy) - 2.5 * SQR (dy) + 1;
-  v[2] = -1.5 * CUBE (dy) + 2 * SQR (dy) + 0.5 * dy;
-  v[3] = 0.5 * CUBE (dy) - 0.5 * SQR (dy);
-
-  w[0] = -0.5 * CUBE (dz) + SQR (dz) - 0.5 * dz;
-  w[1] = 1.5 * CUBE (dz) - 2.5 * SQR (dz) + 1;
-  w[2] = -1.5 * CUBE (dz) + 2 * SQR (dz) + 0.5 * dz;
-  w[3] = 0.5 * CUBE (dz) - 0.5 * SQR (dz);
-
-  for (k = 0; k < 4; k++)
-  {
-    q[k] = 0;
-    for (j = 0; j < 4; j++)
-    {
-      r[j] = 0;
-      for (i = 0; i < 4; i++)
-      {
-        r[j] += u[i] * *pv;
-        pv++;
-      }
-      q[k] += v[j] * r[j];
-      pv += xDim - 4;
-    }
-    vox += w[k] * q[k];
-    pv += xyDim - 4 * xDim;
-  }
-  return vox;
-}
diff --git a/gradunwarp/core/legendre_ext.c b/gradunwarp/core/legendre_ext.c
deleted file mode 100644
index 372c629..0000000
--- a/gradunwarp/core/legendre_ext.c
+++ /dev/null
@@ -1,169 +0,0 @@
-/*
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-#
-#   See COPYING file distributed along with the gradunwarp package for the
-#   copyright and license terms.
-#
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-*/
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include <math.h>
-
-/* compile command
- gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.7 -o legendre_ext.so legendre_ext.c
- */
-static PyObject *legendre(PyObject *self, PyObject *args);
-long odd_factorial(int k);
-
-// what function are exported
-static PyMethodDef legendremethods[] = {
-	{"_legendre", legendre, METH_VARARGS},
-	{NULL, NULL}
-};
-
-// This function is essential for an extension for Numpy created in C
-#if PY_MAJOR_VERSION >= 3
-
-static struct PyModuleDef moduledef = {
-        PyModuleDef_HEAD_INIT,
-        "legendre_ext",
-        NULL,
-        -1,
-        legendremethods,
-        NULL,
-        NULL,
-        NULL,
-        NULL
-};
-
-#define INITERROR return NULL
-
-PyObject *
-PyInit_legendre_ext(void)
-
-#else
-#define INITERROR return
-
-void
-initlegendre_ext(void)
-#endif
-{
-#if PY_MAJOR_VERSION >= 3
-    PyObject *module = PyModule_Create(&moduledef);
-#else
-    PyObject *module = Py_InitModule("legendre_ext", legendremethods);
-#endif
-
-    if (module == NULL)
-        INITERROR;
-
-    import_array();
-
-#if PY_MAJOR_VERSION >= 3
-    return module;
-#endif
-}
-
-// the data should be FLOAT32 and should be ensured in the wrapper 
-static PyObject *legendre(PyObject *self, PyObject *args)
-{
-	PyArrayObject *x, *result;
-    long nu, mu;
-
-	// We expect 1 argument of the PyArray_Type
-    // and 2 of Python float32
-	if(!PyArg_ParseTuple(args, "llO!", 
-                &nu, &mu,
-				&PyArray_Type, &x)) return NULL;
-
-	if ( NULL == x ) return NULL;
-
-    // Error check
-    if ( mu < 0.0 || mu > nu )
-    {
-        printf("Error!: require legendre computation to have 0 <= mu <=nu,\n");
-        printf("but mu=%d and nu=%d\n", mu, nu);
-        return NULL; 
-    }
-
-	//result matrix is the same size as x and is float
-	result = (PyArrayObject*) PyArray_ZEROS(PyArray_NDIM(x), x->dimensions, NPY_FLOAT, 0); 
-	// This is for reference counting ( I think )
-	PyArray_FLAGS(result) |= NPY_OWNDATA; 
-
-    PyArrayIterObject *itr_x, *itr_r;
-    int s, n;
-    float *px, *pr, p_nu, p_nu_prev;
-
-    itr_x = (PyArrayIterObject *) PyArray_IterNew(x);
-    itr_r = (PyArrayIterObject *) PyArray_IterNew(result);
-    while(PyArray_ITER_NOTDONE(itr_x)) 
-    {
-		px = (float *) PyArray_ITER_DATA(itr_x);
-        // unfortunate error check in the main loop
-        if ( abs(*px) > 1.0 )
-        {
-            printf("Error! require -1 <= x <= 1 in the legendre computation.\n");
-            printf("but got x=%f\n", x);
-            return NULL; 
-        }
-		pr = (float *) PyArray_ITER_DATA(itr_r);
-
-        // Compute the initial term in the recursion
-        if ( mu )
-        {
-           s = 1;
-           if ( mu & 1)
-              s = -1;
-              p_nu = s * odd_factorial(2 * mu - 1) * pow(sqrt( 1.0 - *px * *px), mu);
-        }
-        else 
-           p_nu = 1.0;
-
-        // special case.. clear up and return
-        if ( mu == nu )
-        {
-            *pr = p_nu;
-		    PyArray_ITER_NEXT(itr_x);
-		    PyArray_ITER_NEXT(itr_r);
-            continue;
-        }
-
-        // compute the next term in recursion
-        p_nu_prev = p_nu;
-        p_nu = *px * (2 * mu + 1) * p_nu;
-
-        // special case.. clear up and return
-        if ( nu == mu + 1)
-        {
-            *pr = p_nu;
-		    PyArray_ITER_NEXT(itr_x);
-		    PyArray_ITER_NEXT(itr_r);
-            continue;
-        }
-
-        for(n=mu+2; n<nu+1; n++)
-        {
-            *pr = (*px * (2 * n - 1) * p_nu - (n + mu - 1) * p_nu_prev) / (n - mu);
-            p_nu_prev = p_nu;
-            p_nu = *pr;
-        }
-
-		PyArray_ITER_NEXT(itr_x);
-		PyArray_ITER_NEXT(itr_r);
-    }
-
-    return result;
-}
-
-long odd_factorial(int k)
-{
-    long f = (long)k;
-    while ( k >= 3)
-    {
-        k -= 2;
-        f *= k;
-    }
-    return f;
-}
diff --git a/gradunwarp/core/tests/test_utils.py b/gradunwarp/core/tests/test_utils.py
deleted file mode 100644
index 46161b1..0000000
--- a/gradunwarp/core/tests/test_utils.py
+++ /dev/null
@@ -1,46 +0,0 @@
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-#
-#   See COPYING file distributed along with the NiBabel package for the
-#   copyright and license terms.
-#
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-'''Test multiple utils'''
-import os
-import numpy as np
-from numpy.testing import assert_equal, assert_array_equal, \
-    assert_array_almost_equal, assert_almost_equal
-from gradunwarp.core.utils import interp3
-from gradunwarp.core.utils import legendre
-
-
-def test_interp3():
-    arr = np.linspace(-4, 4, 6000)
-    arr = np.sin(arr)
-    arr = arr.reshape(10, 20, 30).astype('float32')
-
-    # sanity check
-    ex1 = arr[4, 11, 15]
-    ac1 = interp3(arr, np.array([4.]), np.array([11.]), np.array([15.]))
-    assert_almost_equal(ex1, ac1[0], 5)
-
-    ex2 = arr[5, 12, 16]
-    ac2 = interp3(arr, np.array([5.]), np.array([12.]), np.array([16.]))
-    assert_almost_equal(ex2, ac2[0], 5)
-
-    ex3 = np.array([-0.33291185, -0.24946867, -0.1595035,
-                    -0.06506295,  0.03180848, 0.12906794,
-                    0.22467692,  0.31659847,  0.40280011,
-                    0.48125309])
-    gridn = 10
-    R1 = np.linspace(4., 5., gridn).astype('float32')
-    C1 = np.linspace(11., 12., gridn).astype('float32')
-    S1 = np.linspace(15., 16., gridn).astype('float32')
-    ac3 = interp3(arr, R1, C1, S1)
-    assert_array_almost_equal(ex3, ac3)
-
-
-def test_legendre():
-    arr = np.array([0.1, 0.2, 0.3])
-    ex1 = np.array([44.83644, 75.85031, 82.44417])
-    ac1 = legendre(6, 3, arr)
-    assert_array_almost_equal(ex1, ac1, 5)
diff --git a/gradunwarp/core/transform_coordinates_ext.c b/gradunwarp/core/transform_coordinates_ext.c
deleted file mode 100644
index 8323f0b..0000000
--- a/gradunwarp/core/transform_coordinates_ext.c
+++ /dev/null
@@ -1,153 +0,0 @@
-/*
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-#
-#   See COPYING file distributed along with the gradunwarp package for the
-#   copyright and license terms.
-#
-### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ##
-*/
-#include "Python.h"
-#include "numpy/arrayobject.h"
-#include <math.h>
-
-/* compile command
- gcc -shared -pthread -fPIC -fwrapv -O2 -Wall -fno-strict-aliasing -I/usr/include/python2.7 -o transform_coordinates_ext.so transform_coordinates_ext.c
- */
-static PyObject *transform_coordinates(PyObject *self, PyObject *args);
-
-
-// what function are exported
-static PyMethodDef transform_coordinatesmethods[] = {
-	{"_transform_coordinates", transform_coordinates, METH_VARARGS},
-	{NULL, NULL}
-};
-
-
-// This function is essential for an extension for Numpy created in C
-#if PY_MAJOR_VERSION >= 3
-
-static struct PyModuleDef moduledef = {
-        PyModuleDef_HEAD_INIT,
-        "transform_coordinates_ext",
-        NULL,
-        -1,
-        transform_coordinatesmethods,
-        NULL,
-        NULL,
-        NULL,
-        NULL
-};
-
-#define INITERROR return NULL
-
-PyObject *
-PyInit_transform_coordinates_ext(void)
-
-#else
-#define INITERROR return
-
-void
-inittransform_coordinates_ext(void)
-#endif
-{
-#if PY_MAJOR_VERSION >= 3
-    PyObject *module = PyModule_Create(&moduledef);
-#else
-    PyObject *module = Py_InitModule("transform_coordinates_ext", transform_coordinatesmethods);
-#endif
-
-    if (module == NULL)
-        INITERROR;
-
-    import_array();
-
-#if PY_MAJOR_VERSION >= 3
-    return module;
-#endif
-}
-
-
-// the data should be FLOAT32 and should be ensured in the wrapper 
-static PyObject *transform_coordinates(PyObject *self, PyObject *args)
-{
-    float *m;
-    float *x, *y, *z, *xm, *ym, *zm;
-    PyArrayObject *X, *Y, *Z, *mat;
-    PyArrayObject *Xm, *Ym, *Zm;
-	// We expect 4 arguments of the PyArray_Type
-	if(!PyArg_ParseTuple(args, "O!O!O!O!", 
-                &PyArray_Type, &X,
-                &PyArray_Type, &Y,
-                &PyArray_Type, &Z,
-                &PyArray_Type, &mat)) return NULL;
-
-    if ( NULL == X ) return NULL;
-    if ( NULL == Y ) return NULL;
-    if ( NULL == Z ) return NULL;
-    if ( NULL == mat ) return NULL;
-
-    // result matrices are the same size and float
-	Xm = (PyArrayObject*) PyArray_ZEROS(PyArray_NDIM(X), X->dimensions, NPY_FLOAT, 0); 
-	Ym = (PyArrayObject*) PyArray_ZEROS(PyArray_NDIM(X), X->dimensions, NPY_FLOAT, 0); 
-	Zm = (PyArrayObject*) PyArray_ZEROS(PyArray_NDIM(X), X->dimensions, NPY_FLOAT, 0); 
-  
-	// This is for reference counting ( I think )
-	PyArray_FLAGS(Xm) |= NPY_OWNDATA; 
-	PyArray_FLAGS(Ym) |= NPY_OWNDATA; 
-	PyArray_FLAGS(Zm) |= NPY_OWNDATA; 
-
-	// massive use of iterators to progress through the data
-	PyArrayIterObject *itr_x, *itr_y, *itr_z;
-	PyArrayIterObject *itr_xm, *itr_ym, *itr_zm;
-    itr_x = (PyArrayIterObject *) PyArray_IterNew(X);
-    itr_y = (PyArrayIterObject *) PyArray_IterNew(Y);
-    itr_z = (PyArrayIterObject *) PyArray_IterNew(Z);
-    itr_xm = (PyArrayIterObject *) PyArray_IterNew(Xm);
-    itr_ym = (PyArrayIterObject *) PyArray_IterNew(Ym);
-    itr_zm = (PyArrayIterObject *) PyArray_IterNew(Zm);
-    /*m = (float *)PyArray_DATA(mat);
-    printf("%f %f %f %f\n", m[0], m[1], m[2], m[3]);
-    printf("%f %f %f %f\n", m[4], m[5], m[6], m[7]);
-    printf("%f %f %f %f\n", m[8], m[9], m[10], m[11]);
-    */
-    float *m00, *m01, *m02, *m03;
-    float *m10, *m11, *m12, *m13;
-    float *m20, *m21, *m22, *m23;
-    m00 = (float *)PyArray_GETPTR2(mat, 0, 0);
-    m01 = (float *)PyArray_GETPTR2(mat, 0, 1);
-    m02 = (float *)PyArray_GETPTR2(mat, 0, 2);
-    m03 = (float *)PyArray_GETPTR2(mat, 0, 3);
-    m10 = (float *)PyArray_GETPTR2(mat, 1, 0);
-    m11 = (float *)PyArray_GETPTR2(mat, 1, 1);
-    m12 = (float *)PyArray_GETPTR2(mat, 1, 2);
-    m13 = (float *)PyArray_GETPTR2(mat, 1, 3);
-    m20 = (float *)PyArray_GETPTR2(mat, 2, 0);
-    m21 = (float *)PyArray_GETPTR2(mat, 2, 1);
-    m22 = (float *)PyArray_GETPTR2(mat, 2, 2);
-    m23 = (float *)PyArray_GETPTR2(mat, 2, 3);
-
-    // start the iteration
-    while(PyArray_ITER_NOTDONE(itr_x))
-    {
-        x = (float *) PyArray_ITER_DATA(itr_x);
-        y = (float *) PyArray_ITER_DATA(itr_y);
-        z = (float *) PyArray_ITER_DATA(itr_z);
-        xm = (float *) PyArray_ITER_DATA(itr_xm);
-        ym = (float *) PyArray_ITER_DATA(itr_ym);
-        zm = (float *) PyArray_ITER_DATA(itr_zm);
-
-        // transform coordinates
-        *xm = *x * *m00 + *y * *m01 + *z * *m02 + *m03;
-        *ym = *x * *m10 + *y * *m11 + *z * *m12 + *m13;
-        *zm = *x * *m20 + *y * *m21 + *z * *m22 + *m23;
-
-		PyArray_ITER_NEXT(itr_x);
-		PyArray_ITER_NEXT(itr_y);
-		PyArray_ITER_NEXT(itr_z);
-		PyArray_ITER_NEXT(itr_xm);
-		PyArray_ITER_NEXT(itr_ym);
-		PyArray_ITER_NEXT(itr_zm);
-    }
-
- return Py_BuildValue("OOO", Xm, Ym, Zm);
-}
diff --git a/gradunwarp/core/unwarp_resample.py b/gradunwarp/core/unwarp_resample.py
index 89ab54f..46f3bc6 100644
--- a/gradunwarp/core/unwarp_resample.py
+++ b/gradunwarp/core/unwarp_resample.py
@@ -19,23 +19,23 @@
 from .globals import siemens_max_det
 import nibabel as nib
 import subprocess
+import scipy.special
 
-#np.seterr(all='raise')
+# np.seterr(all='raise')
 
-log = logging.getLogger('gradunwarp')
+log = logging.getLogger("gradunwarp")
 
 
 class Unwarper(object):
-    '''
-    '''
-    def __init__(self, vol, m_rcs2ras, vendor, coeffs, fileName):
-        '''
-        '''
-        self.vol = vol
-        self.m_rcs2ras = m_rcs2ras
+    """ """
+
+    def __init__(self, input_nii, vendor, coeffs, fileName):
+        """ """
+        self.input_nii = input_nii
+        self.vol, self.m_rcs2ras = self.input_nii.get_fdata(), self.input_nii.affine
         self.vendor = vendor
         self.coeffs = coeffs
-        self.name=fileName
+        self.name = fileName
         self.warp = False
         self.nojac = False
         self.m_rcs2lai = None
@@ -49,10 +49,10 @@ def __init__(self, vol, m_rcs2ras, vendor, coeffs, fileName):
         self.order = 1
 
     def eval_spharm_grid(self, vendor, coeffs):
-        ''' 
+        """
         We evaluate the spherical harmonics on a less sampled grid.
         This is a spacetime vs accuracy tradeoff.
-        '''
+        """
         # init the grid first
         if not self.fovmin:
             fovmin = globals.siemens_fovmin
@@ -68,69 +68,68 @@ def eval_spharm_grid(self, vendor, coeffs):
             numpoints = self.numpoints
 
         # convert to mm
-        fovmin = fovmin * 1000.
-        fovmax = fovmax * 1000.
+        fovmin = fovmin * 1000.0
+        fovmax = fovmax * 1000.0
         # the grid in meters. this is needed for spherical harmonics
         vec = np.linspace(fovmin, fovmax, numpoints)
         gvx, gvy, gvz = utils.meshgrid(vec, vec, vec)
         # mm
         cf = (fovmax - fovmin) / numpoints
-        
+
         # deduce the transformation from rcs to grid
-        g_rcs2xyz = np.array( [[0, cf, 0, fovmin],
-                               [cf, 0, 0, fovmin],
-                               [0, 0, cf, fovmin],
-                               [0, 0, 0, 1]], dtype=np.float32 )
+        g_rcs2xyz = np.array(
+            [[0, cf, 0, fovmin], [cf, 0, 0, fovmin], [0, 0, cf, fovmin], [0, 0, 0, 1]],
+            dtype=np.float32,
+        )
 
         # get the grid to rcs transformation also
         g_xyz2rcs = np.linalg.inv(g_rcs2xyz)
 
         # indices into the gradient displacement vol
-        gr, gc, gs = utils.meshgrid(np.arange(numpoints), np.arange(numpoints),
-                                 np.arange(numpoints), dtype=np.float32)
-
-        log.info('Evaluating spherical harmonics')
-        log.info('on a ' + str(numpoints) + '^3 grid')
-        log.info('with extents ' + str(fovmin) + 'mm to ' + str(fovmax) + 'mm')
+        gr, gc, gs = utils.meshgrid(
+            np.arange(numpoints),
+            np.arange(numpoints),
+            np.arange(numpoints),
+            dtype=np.float32,
+        )
+
+        log.info("Evaluating spherical harmonics")
+        log.info("on a " + str(numpoints) + "^3 grid")
+        log.info("with extents " + str(fovmin) + "mm to " + str(fovmax) + "mm")
         gvxyz = CV(gvx, gvy, gvz)
         _dv, _dxyz = eval_spherical_harmonics(coeffs, vendor, gvxyz)
-            
-        return CV(_dv.x, _dv.y, _dv.z), g_xyz2rcs
 
+        return CV(_dv.x, _dv.y, _dv.z), g_xyz2rcs
 
     def run(self):
-        '''
-        '''
-        #pdb.set_trace()
+        """ """
+        # pdb.set_trace()
         # define polarity based on the warp requested
-        self.polarity = 1.
+        self.polarity = 1.0
         if self.warp:
-            self.polarity = -1.
+            self.polarity = -1.0
 
-        # Evaluate spherical harmonics on a smaller grid 
+        # Evaluate spherical harmonics on a smaller grid
         dv, g_xyz2rcs = self.eval_spharm_grid(self.vendor, self.coeffs)
 
         # transform RAS-coordinates into LAI-coordinates
-        m_ras2lai = np.array([[-1.0, 0.0, 0.0, 0.0],
-                             [0.0, 1.0, 0.0, 0.0],
-                             [0.0, 0.0, -1.0, 0.0],
-                             [0.0, 0.0, 0.0, 1.0]], dtype=np.float64)
+        m_ras2lai = np.array(
+            [
+                [-1.0, 0.0, 0.0, 0.0],
+                [0.0, 1.0, 0.0, 0.0],
+                [0.0, 0.0, -1.0, 0.0],
+                [0.0, 0.0, 0.0, 1.0],
+            ],
+            dtype=np.float64,
+        )
         m_rcs2lai = np.dot(m_ras2lai, self.m_rcs2ras)
         m_rcs2lai_nohalf = m_rcs2lai[:, :]
 
-        # indices of image volume
-        '''
-        nr, nc, ns = self.vol.shape[:3]
-        vc3, vr3, vs3 = utils.meshgrid(np.arange(nr), np.arange(nc), np.arange(ns), dtype=np.float32)
-        vrcs = CV(x=vr3, y=vc3, z=vs3)
-        vxyz = utils.transform_coordinates(vrcs, m_rcs2lai)
-        '''
-
         # account for half-voxel shift in R and C directions
         halfvox = np.zeros((4, 4))
         halfvox[0, 3] = m_rcs2lai[0, 0] / 2.0
         halfvox[1, 3] = m_rcs2lai[1, 1] / 2.0
-        #m_rcs2lai = m_rcs2lai + halfvox
+        # m_rcs2lai = m_rcs2lai + halfvox
 
         # extract rotational and scaling parts of the transformation matrix
         # ignore the translation part
@@ -145,16 +144,20 @@ def run(self):
         # The differentials can be determined by mapping a vector of 1s
         # through rotation and scaling, where any mirror
         # will impose a negation
-        ones = CV(1., 1., 1.)
-        dxyz = utils.transform_coordinates_old(ones, r_rcs2lai)
+        ones = CV(1.0, 1.0, 1.0)
+        dxyz = utils.transform_coordinates(ones, r_rcs2lai)
 
         # do the nonlinear unwarp
-        if self.vendor == 'siemens':
-            self.out, self.vjacout = self.non_linear_unwarp_siemens(self.vol.shape, dv, dxyz,
-                                                                 m_rcs2lai, m_rcs2lai_nohalf, g_xyz2rcs)
-
-    def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_nohalf, g_xyz2rcs):
-        ''' Performs the crux of the unwarping.
+        if self.vendor == "siemens":
+            self.out, self.vjacout = self.non_linear_unwarp_siemens(
+                self.vol.shape, dv, dxyz, m_rcs2lai, m_rcs2lai_nohalf, g_xyz2rcs
+            )
+        del dv, g_xyz2rcs, self.vjacout
+
+    def non_linear_unwarp_siemens(
+        self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_nohalf, g_xyz2rcs
+    ):
+        """Performs the crux of the unwarping.
         It's agnostic to Siemens or GE and uses more functions to
         do the processing separately.
 
@@ -174,8 +177,8 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh
             x,y and z coordinates of the unwarped coordinates
         vjacmult_lps : np.array
             the jacobian multiplier (determinant)
-        '''
-        log.info('Evaluating the jacobian multiplier')
+        """
+        log.info("Evaluating the jacobian multiplier")
         nr, nc, ns = self.vol.shape[:3]
         if not self.nojac:
             jim2 = np.zeros((nr, nc), dtype=np.float32)
@@ -185,7 +188,7 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh
             else:
                 vjacdet_lps = eval_siemens_jacobian_mult(dv, dxyz)
 
-        # essentially pre-allocating everything 
+        # essentially pre-allocating everything
         out = np.zeros((nr, nc, ns), dtype=np.float32)
         fullWarp = np.zeros((nr, nc, ns, 3), dtype=np.float32)
 
@@ -199,159 +202,164 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh
         vc, vr = utils.meshgrid(np.arange(nc), np.arange(nr))
 
         # Compute transform to map the internal voxel coordinates to FSL scaled mm coordinates
-        try:
-            pixdim1=float((subprocess.Popen(['fslval', self.name,'pixdim1'], stdout=subprocess.PIPE).communicate()[0]).strip())
-        except ValueError:
-            log.error('Failure during fslval call. Make sure fslval, fslhd, fslorient are in your PATH, and that FSLOUTPUTTYPE is set.')
-            sys.exit(1)
-        
-        pixdim2=float((subprocess.Popen(['fslval', self.name,'pixdim2'], stdout=subprocess.PIPE).communicate()[0]).strip())
-        pixdim3=float((subprocess.Popen(['fslval', self.name,'pixdim3'], stdout=subprocess.PIPE).communicate()[0]).strip())
-        dim1=float((subprocess.Popen(['fslval', self.name,'dim1'], stdout=subprocess.PIPE).communicate()[0]).strip())
-        outputOrient=subprocess.Popen(['fslorient', self.name], stdout=subprocess.PIPE).communicate()[0].strip()
-        if outputOrient == b'NEUROLOGICAL':
-            log.info('Input volume is NEUROLOGICAL orientation. Flipping x-axis in output fullWarp_abs.nii.gz')
+        pixdim1, pixdim2, pixdim3 = self.input_nii.header.get_zooms()[:3]
+        dim1 = self.input_nii.shape[0]
+
+        # TODO: figure out if that logic works
+        if np.linalg.det(self.input_nii.affine) > 0:
+            log.info(
+                "Input volume is NEUROLOGICAL orientation. Flipping x-axis in output fullWarp_abs.nii.gz"
+            )
             # if neurological then flip x coordinate (both here in premat and later in postmat)
-            m_vox2fsl = np.array([[-1.0*pixdim1, 0.0, 0.0, pixdim1*(dim1-1)],
-                              [0.0, pixdim2, 0.0, 0.0],
-                              [0.0, 0.0, pixdim3, 0.0],
-                              [0.0, 0.0, 0.0, 1.0]], dtype=np.float64)
+            m_vox2fsl = np.array(
+                [
+                    [-1.0 * pixdim1, 0.0, 0.0, pixdim1 * (dim1 - 1)],
+                    [0.0, pixdim2, 0.0, 0.0],
+                    [0.0, 0.0, pixdim3, 0.0],
+                    [0.0, 0.0, 0.0, 1.0],
+                ],
+                dtype=np.float64,
+            )
         else:
-            m_vox2fsl = np.array([[pixdim1, 0.0, 0.0, 0.0],
-                              [0.0, pixdim2, 0.0, 0.0],
-                              [0.0, 0.0, pixdim3, 0.0],
-                              [0.0, 0.0, 0.0, 1.0]], dtype=np.float64)
-															
-        log.info('Unwarping slice by slice')
+            m_vox2fsl = np.array(
+                [
+                    [pixdim1, 0.0, 0.0, 0.0],
+                    [0.0, pixdim2, 0.0, 0.0],
+                    [0.0, 0.0, pixdim3, 0.0],
+                    [0.0, 0.0, 0.0, 1.0],
+                ],
+                dtype=np.float64,
+            )
+
+        log.info("Unwarping slice by slice")
         # for every slice
         for s in range(ns):
             # pretty print
             sys.stdout.flush()
-            if (s+1) % 10 == 0:
-                print(s+1, end=' ')
+            if (s + 1) % 10 == 0:
+                print(s + 1, end=" ")
             else:
-                print('.', end=' ')
-                
+                print(".", end=" ")
+
             # hopefully, free memory
             gc.collect()
             # init to 0
-            dvx.fill(0.)
-            dvy.fill(0.)
-            dvz.fill(0.)
-            im_.fill(0.)
-            
+            dvx.fill(0.0)
+            dvy.fill(0.0)
+            dvz.fill(0.0)
+            im_.fill(0.0)
+
             vs = np.ones(vr.shape) * s
             vrcs = CV(vr, vc, vs)
             vxyz = utils.transform_coordinates(vrcs, m_rcs2lai_nohalf)
             vrcsg = utils.transform_coordinates(vxyz, g_xyz2rcs)
-            ndimage.interpolation.map_coordinates(dv.x,
-                                                  vrcsg,
-                                                  output=dvx,
-                                                  order=self.order)
-            ndimage.interpolation.map_coordinates(dv.y,
-                                                  vrcsg,
-                                                  output=dvy,
-                                                  order=self.order)
-            ndimage.interpolation.map_coordinates(dv.z,
-                                                  vrcsg,
-                                                  output=dvz,
-                                                  order=self.order)
+            ndimage.interpolation.map_coordinates(
+                dv.x, vrcsg, output=dvx, order=self.order
+            )
+            ndimage.interpolation.map_coordinates(
+                dv.y, vrcsg, output=dvy, order=self.order
+            )
+            ndimage.interpolation.map_coordinates(
+                dv.z, vrcsg, output=dvz, order=self.order
+            )
             # new locations of the image voxels in XYZ ( LAI ) coords
 
-            #dvx.fill(0.)
-            #dvy.fill(0.)
-            #dvz.fill(0.)
-            
-            vxyzw = CV(x=vxyz.x + self.polarity * dvx,
-                       y=vxyz.y + self.polarity * dvy,
-                       z=vxyz.z + self.polarity * dvz)
+            # dvx.fill(0.)
+            # dvy.fill(0.)
+            # dvz.fill(0.)
+
+            vxyzw = CV(
+                x=vxyz.x + self.polarity * dvx,
+                y=vxyz.y + self.polarity * dvy,
+                z=vxyz.z + self.polarity * dvz,
+            )
 
             # convert the locations got into RCS indices
-            vrcsw = utils.transform_coordinates(vxyzw,
-                                                np.linalg.inv(m_rcs2lai))
+            vrcsw = utils.transform_coordinates(vxyzw, np.linalg.inv(m_rcs2lai))
             # map the internal voxel coordinates to FSL scaled mm coordinates
             vfsl = utils.transform_coordinates(vrcsw, m_vox2fsl)
 
-
-            #im_ = utils.interp3(self.vol, vrcsw.x, vrcsw.y, vrcsw.z)
-            ndimage.interpolation.map_coordinates(self.vol,
-                                                  vrcsw,
-                                                  output=im_,
-                                                  order=self.order)
+            # im_ = utils.interp3(self.vol, vrcsw.x, vrcsw.y, vrcsw.z)
+            ndimage.interpolation.map_coordinates(
+                self.vol, vrcsw, output=im_, order=self.order
+            )
             # find NaN voxels, and set them to 0
-            im_[np.where(np.isnan(im_))] = 0.
-            im_[np.where(np.isinf(im_))] = 0.
+            im_[np.where(np.isnan(im_))] = 0.0
+            im_[np.where(np.isinf(im_))] = 0.0
             im2[vr, vc] = im_
 
-            #img = nib.Nifti1Image(dvx,np.eye(4))
-            #nib.save(img,"x"+str(s).zfill(3)+".nii.gz")
-            #img = nib.Nifti1Image(dvy,np.eye(4))
-            #nib.save(img,"y"+str(s).zfill(3)+".nii.gz")
-            #img = nib.Nifti1Image(dvz,np.eye(4))
-            #nib.save(img,"z"+str(s).zfill(3)+".nii.gz")
+            # img = nib.Nifti1Image(dvx,np.eye(4))
+            # nib.save(img,"x"+str(s).zfill(3)+".nii.gz")
+            # img = nib.Nifti1Image(dvy,np.eye(4))
+            # nib.save(img,"y"+str(s).zfill(3)+".nii.gz")
+            # img = nib.Nifti1Image(dvz,np.eye(4))
+            # nib.save(img,"z"+str(s).zfill(3)+".nii.gz")
 
             # Multiply the intensity with the Jacobian det, if needed
             if not self.nojac:
-                vjacdet_lpsw.fill(0.)
-                jim2.fill(0.)
+                vjacdet_lpsw.fill(0.0)
+                jim2.fill(0.0)
                 # if polarity is negative, the jacobian is also inversed
                 if self.polarity == -1:
-                    vjacdet_lps = 1. / vjacdet_lps
-
-                ndimage.interpolation.map_coordinates(vjacdet_lps,
-                                                      vrcsg,
-                                                      output=vjacdet_lpsw,
-                                                      order=self.order)
-                vjacdet_lpsw[np.where(np.isnan(vjacdet_lpsw))] = 0.
-                vjacdet_lpsw[np.where(np.isinf(vjacdet_lpsw))] = 0.
+                    vjacdet_lps = 1.0 / vjacdet_lps
+
+                ndimage.interpolation.map_coordinates(
+                    vjacdet_lps, vrcsg, output=vjacdet_lpsw, order=self.order
+                )
+                vjacdet_lpsw[np.where(np.isnan(vjacdet_lpsw))] = 0.0
+                vjacdet_lpsw[np.where(np.isinf(vjacdet_lpsw))] = 0.0
                 jim2[vr, vc] = vjacdet_lpsw
                 im2 = im2 * jim2
                 vjacout[..., s] = jim2
 
-            fullWarp[...,s,0]=vfsl.x
-            fullWarp[...,s,1]=vfsl.y
-            fullWarp[...,s,2]=vfsl.z
+            fullWarp[..., s, 0] = vfsl.x
+            fullWarp[..., s, 1] = vfsl.y
+            fullWarp[..., s, 2] = vfsl.z
             out[..., s] = im2
 
         print()
-       
-        img=nib.Nifti1Image(fullWarp,self.m_rcs2ras)
-        nib.save(img,"fullWarp_abs.nii.gz")
+
+        self.fullWarp = nib.Nifti1Image(fullWarp, self.m_rcs2ras)
         # return image and the jacobian
+        del vrcsw, vfsl, vxyzw, vrcs, vxyz, vrcsg, fullWarp
         return out, vjacout
 
     def write(self, outfile):
-        log.info('Writing output to ' + outfile)
+        log.info("Writing output to " + outfile)
         # if out datatype is float64 make it float32
         if self.out.dtype == np.float64:
             self.out = self.out.astype(np.float32)
-        if outfile.endswith('.nii') or outfile.endswith('.nii.gz'):
+        if outfile.endswith(".nii") or outfile.endswith(".nii.gz"):
             img = nib.Nifti1Image(self.out, self.m_rcs2ras)
-        if outfile.endswith('.mgh') or outfile.endswith('.mgz'):
-            #self.out = self.out.astype(self.vol.dtype)
+        if outfile.endswith(".mgh") or outfile.endswith(".mgz"):
+            # self.out = self.out.astype(self.vol.dtype)
             img = nib.MGHImage(self.out, self.m_rcs2ras)
         nib.save(img, outfile)
 
+        suffix = outfile.split('_')[-1]
+        nib.save(self.fullWarp, outfile.replace(suffix, f"desc-{suffix.replace('.nii.gz','')}_warp.nii.gz"))
+
 
 def eval_siemens_jacobian_mult(F, dxyz):
-    '''
-    '''
+    """ """
     d0, d1, d2 = dxyz.x, dxyz.y, dxyz.z
-    #print F.x.shape, d0, d1, d2
+    # print F.x.shape, d0, d1, d2
 
     if d0 == 0 or d1 == 0 or d2 == 0:
-        raise ValueError('weirdness found in Jacobian calculation')
+        raise ValueError("weirdness found in Jacobian calculation")
 
     dFxdx, dFxdy, dFxdz = np.gradient(F.x, d0, d1, d2)
     dFydx, dFydy, dFydz = np.gradient(F.y, d0, d1, d2)
     dFzdx, dFzdy, dFzdz = np.gradient(F.z, d0, d1, d2)
 
-    jacdet = (1. + dFxdx) * (1. + dFydy) * (1. + dFzdz) \
-           - (1. + dFxdx) * dFydz * dFzdy \
-           - dFxdy * dFydx * (1. + dFzdz) \
-           + dFxdy * dFydz * dFzdx \
-           + dFxdz * dFydx * dFzdy \
-           - dFxdz * (1. + dFydy) * dFzdx
+    jacdet = (
+        (1.0 + dFxdx) * (1.0 + dFydy) * (1.0 + dFzdz)
+        - (1.0 + dFxdx) * dFydz * dFzdy
+        - dFxdy * dFydx * (1.0 + dFzdz)
+        + dFxdy * dFydz * dFzdx
+        + dFxdz * dFydx * dFzdy
+        - dFxdz * (1.0 + dFydy) * dFzdx
+    )
     jacdet = np.abs(jacdet)
     jacdet[np.where(jacdet > siemens_max_det)] = siemens_max_det
 
@@ -359,7 +367,7 @@ def eval_siemens_jacobian_mult(F, dxyz):
 
 
 def eval_spherical_harmonics(coeffs, vendor, vxyz):
-    ''' Evaluate spherical harmonics
+    """Evaluate spherical harmonics
 
     Parameters
     ----------
@@ -371,39 +379,38 @@ def eval_spherical_harmonics(coeffs, vendor, vxyz):
         in the function
     resolution : float
         (optional) useful in case vxyz is scalar
-    '''
+    """
     # convert radius into mm
-    R0 = coeffs.R0_m  * 1000
+    R0 = coeffs.R0_m * 1000
 
     x, y, z = vxyz
 
-    #pdb.set_trace()
+    # pdb.set_trace()
     # log.info('calculating displacements (mm) '
     #        'using spherical harmonics coeffcients...')
-    if vendor == 'siemens':
-        log.info('along x...')
+    if vendor == "siemens":
+        log.info("along x...")
         bx = siemens_B(coeffs.alpha_x, coeffs.beta_x, x, y, z, R0)
-        log.info('along y...')
+        log.info("along y...")
         by = siemens_B(coeffs.alpha_y, coeffs.beta_y, x, y, z, R0)
-        log.info('along z...')
+        log.info("along z...")
         bz = siemens_B(coeffs.alpha_z, coeffs.beta_z, x, y, z, R0)
     else:
         # GE
-        log.info('along x...')
+        log.info("along x...")
         bx = ge_D(coeffs.alpha_x, coeffs.beta_x, x, y, z)
-        log.info('along y...')
+        log.info("along y...")
         by = ge_D(coeffs.alpha_y, coeffs.beta_y, x, y, z)
-        log.info('along z...')
+        log.info("along z...")
         bz = siemens_B(coeffs.alpha_z, coeffs.beta_z, x, y, z, R0)
         bz = ge_D(coeffs.alpha_z, coeffs.beta_z, x, y, z)
 
     return CV(bx * R0, by * R0, bz * R0), CV(x, y, z)
 
 
-#@profile
+# @profile
 def siemens_B(alpha, beta, x1, y1, z1, R0):
-    ''' Calculate displacement field from Siemens coefficients
-    '''
+    """Calculate displacement field from Siemens coefficients"""
     nmax = alpha.shape[0] - 1
     x1 = x1 + 0.0001  # hack to avoid singularities at R=0
 
@@ -417,22 +424,22 @@ def siemens_B(alpha, beta, x1, y1, z1, R0):
         f = np.power(r / R0, n)
         for m in range(0, n + 1):
             f2 = alpha[n, m] * np.cos(m * phi) + beta[n, m] * np.sin(m * phi)
-            _ptemp = utils.legendre(n, m, np.cos(theta))
-            #_ptemp = scipy.special.lpmv(m, n, np.cos(theta))
+            # _ptemp = utils.legendre(n, m, np.cos(theta))
+            _p = scipy.special.lpmv(m, n, np.cos(theta))
             normfact = 1
             # this is Siemens normalization
             if m > 0:
-                normfact = math.pow(-1, m) * \
-                math.sqrt(float((2 * n + 1) * factorial(n - m)) \
-                          / float(2 * factorial(n + m)))
-            _p = normfact * _ptemp
-            b = b + f * _p * f2
+                normfact = math.pow(-1, m) * math.sqrt(
+                    float((2 * n + 1) * factorial(n - m)) / float(2 * factorial(n + m))
+                )
+            _p *= normfact
+            b += f * _p * f2
     return b
 
 
 def ge_D(alpha, beta, x1, y1, z1):
-    ''' GE Gradwarp coeffs define the error rather than the total
-    gradient field'''
+    """GE Gradwarp coeffs define the error rather than the total
+    gradient field"""
 
     nmax = alpha.shape[0] - 1
     x1 = x1 + 0.0001  # hack to avoid singularities
@@ -448,9 +455,27 @@ def ge_D(alpha, beta, x1, y1, z1):
         # So GE uses the usual unnormalized legendre polys.
         f = np.power(r, n)
         for m in range(0, n + 1):
-            f2 = alpha[n, m] * np.cos(m * theta) + beta[n, m] \
-            * np.sin(m * theta)
-            _p = utils.legendre(n, m, np.cos(phi))
-            d = d + f * _p * f2
+            f2 = alpha[n, m] * np.cos(m * theta) + beta[n, m] * np.sin(m * theta)
+            _p = scipy.special.lpmv(m, n, np.cos(phi))
+            d += f * _p * f2
     d = d / 100.0  # cm back to meters
     return d
+
+
+
+def get_fslorient(nii):
+    h = nii.header
+    sform = h.get_sform()
+    sform_det = np.linalg.det(sform[:3,:3]) if sform is not None else None
+
+    qform = h.get_qform()
+    qform_det = np.linalg.det(qform[:3,:3]) if qform is not None else None
+
+    det = sform_det or qform_det
+
+    if not sform_det and not qform_det:
+        return 'UNKNOWN'
+    elif sform_det and qform_det and (sform_det * qform_det) < 0 :
+        return 'INCONSISTENT'
+    else:
+        return 'RADIOLOGICAL' if det < 0 else 'NEUROLOGICAL'
diff --git a/gradunwarp/core/utils.py b/gradunwarp/core/utils.py
index 0578352..58905d2 100644
--- a/gradunwarp/core/utils.py
+++ b/gradunwarp/core/utils.py
@@ -16,16 +16,16 @@
 # of a meshgrid belongs to this
 # x, y, z = meshgrid(np.arange(5), np.arange(6), np.arange(7))
 # cv = CoordsVector(x=x, y=y, z=z)
-CoordsVector = namedtuple('CoordsVector', 'x, y, z')
+CoordsVector = namedtuple("CoordsVector", "x, y, z")
 
 
 # this method is deprecated because it's slow and my suspicion that
-# the matrix expressions create unnecessary temp matrices which 
+# the matrix expressions create unnecessary temp matrices which
 # are costly for huge matrices
 def transform_coordinates_old(A, M):
-    ''' 4x4 matrix M operates on orthogonal coordinates arrays
+    """4x4 matrix M operates on orthogonal coordinates arrays
     A1, A2, A3 to give B1, B2, B3
-    '''
+    """
     A1 = A.x
     A2 = A.y
     A3 = A.z
@@ -34,10 +34,11 @@ def transform_coordinates_old(A, M):
     B3 = A1 * M[2, 0] + A2 * M[2, 1] + A3 * M[2, 2] + M[2, 3]
     return CoordsVector(B1, B2, B3)
 
-def transform_coordinates(A, M):
-    ''' 4x4 matrix M operates on orthogonal coordinates arrays
+
+def transform_coordinates_old2(A, M):
+    """4x4 matrix M operates on orthogonal coordinates arrays
     A1, A2, A3 to give B1, B2, B3
-    '''
+    """
     A1 = A.x
     A2 = A.y
     A3 = A.z
@@ -48,21 +49,35 @@ def transform_coordinates(A, M):
     try:
         from .transform_coordinates_ext import _transform_coordinates
     except ImportError:
-        raise ImportError('The transform_coordinates C extension module is missing.' \
-                           ' Fallback code not yet implemented.')
+        raise ImportError(
+            "The transform_coordinates C extension module is missing."
+            " Fallback code not yet implemented."
+        )
 
     B1, B2, B3 = _transform_coordinates(A1, A2, A3, M)
     return CoordsVector(B1, B2, B3)
 
 
+from nibabel.affines import apply_affine
+
+
+def transform_coordinates(A, M):
+    vecs = np.asanyarray([A.x, A.y, A.z]).T
+    vecs_trans = apply_affine(M, vecs)
+    return CoordsVector(
+        vecs_trans[..., 0].T, vecs_trans[..., 1].T, vecs_trans[..., 2].T
+    )
+
+
 def get_vol_affine(infile):
     try:
         import nibabel as nib
     except ImportError:
-        raise ImportError('gradunwarp needs nibabel for I/O of mgz/nifti files.'
-                          ' Please install')
+        raise ImportError(
+            "gradunwarp needs nibabel for I/O of mgz/nifti files." " Please install"
+        )
     nibimage = nib.load(infile)
-    return nibimage.get_data(), nibimage.affine
+    return nibimage.get_fdata(), nibimage.affine
 
 
 # memoized factorial
@@ -76,7 +91,8 @@ def __call__(self, *args):
             self.memo[args] = self.f(*args)
         return self.memo[args]
 
-#factorial = Memoize(math.factorial)
+
+# factorial = Memoize(math.factorial)
 factorial = math.factorial
 
 
@@ -168,28 +184,27 @@ def meshgrid(*xi, **kwargs):
     >>> xx, yy = meshgrid(x, y, sparse=True)
     >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2)
     """
-    copy = kwargs.get('copy', True)
+    copy = kwargs.get("copy", True)
     args = np.atleast_1d(*xi)
     if not isinstance(args, list):
         if args.size > 0:
             return args.copy() if copy else args
         else:
-            raise TypeError('meshgrid() take 1 or more arguments (0 given)')
+            raise TypeError("meshgrid() take 1 or more arguments (0 given)")
 
-    sparse = kwargs.get('sparse', False)
-    indexing = kwargs.get('indexing', 'xy')  # 'ij'
+    sparse = kwargs.get("sparse", False)
+    indexing = kwargs.get("indexing", "xy")  # 'ij'
 
     ndim = len(args)
     s0 = (1,) * ndim
-    output = [x.reshape(s0[:i] + (-1, ) + s0[i + 1::]) \
-              for i, x in enumerate(args)]
+    output = [x.reshape(s0[:i] + (-1,) + s0[i + 1 : :]) for i, x in enumerate(args)]
 
     shape = [x.size for x in output]
 
-    if indexing == 'xy':
+    if indexing == "xy":
         # switch first and second axis
-        output[0].shape = (1, -1) + (1, ) * (ndim - 2)
-        output[1].shape = (-1, 1) + (1, ) * (ndim - 2)
+        output[0].shape = (1, -1) + (1,) * (ndim - 2)
+        output[1].shape = (-1, 1) + (1,) * (ndim - 2)
         shape[0], shape[1] = shape[1], shape[0]
 
     if sparse:
@@ -211,7 +226,7 @@ def ndgrid(*args, **kwargs):
     Same as calling meshgrid with indexing='ij' (see meshgrid for
     documentation).
     """
-    kwargs['indexing'] = 'ij'
+    kwargs["indexing"] = "ij"
     return meshgrid(*args, **kwargs)
 
 
@@ -231,8 +246,7 @@ def legendre_old(nu, mu, x):
     (Abramowitz & Stegun, Section 8.5.)
     """
     if mu < 0 or mu > nu:
-        raise ValueError('require 0 <= mu <= nu, but mu=%d and nu=%d' \
-                         % (nu, mu))
+        raise ValueError("require 0 <= mu <= nu, but mu=%d and nu=%d" % (nu, mu))
     # if abs(x) > 1:
     #    raise ValueError('require -1 <= x <= 1, but x=%f', x)
 
@@ -243,8 +257,8 @@ def legendre_old(nu, mu, x):
         s = 1
         if mu & 1:
             s = -1
-        z = sqrt(1 - x ** 2)
-        p_nu = s * odd_factorial(2 * mu - 1) * z ** mu
+        z = sqrt(1 - x**2)
+        p_nu = s * odd_factorial(2 * mu - 1) * z**mu
 
     if mu == nu:
         return p_nu
@@ -269,24 +283,28 @@ def legendre(nu, mu, x):
     try:
         from .legendre_ext import _legendre
     except ImportError:
-        raise ImportError('The legendre C extension module is missing.' \
-                           ' Fallback legendre code not yet implemented.')
+        raise ImportError(
+            "The legendre C extension module is missing."
+            " Fallback legendre code not yet implemented."
+        )
     nu = int(nu)
     mu = int(mu)
     x = x.astype(np.float32)
     b = _legendre(nu, mu, x)
     return b
-    
+
 
 def interp3(vol, R, C, S):
-    '''
+    """
     TODO
-    '''
+    """
     try:
         from .interp3_ext import _interp3
     except ImportError:
-        raise ImportError('The interp3 C extension module is missing.' \
-                           ' Fallback interp3 code not yet implemented.')
+        raise ImportError(
+            "The interp3 C extension module is missing."
+            " Fallback interp3 code not yet implemented."
+        )
     vol = vol.astype(np.float32)
     vol = np.ascontiguousarray(vol)
     R = R.astype(np.float32)
@@ -299,29 +317,31 @@ def interp3(vol, R, C, S):
     return V
 
 
-if __name__ == '__main__':
-    print('Testing meshgrid...')
-    #import doctest
-    #doctest.testmod()
+if __name__ == "__main__":
+    print("Testing meshgrid...")
+    # import doctest
+    # doctest.testmod()
 
-    print('Profiling interp3...')
-    print('Interpolation for a million coordinates should be' \
-          ' in the order of seconds. Anything significantly less ' \
-          'is considered slow')
+    print("Profiling interp3...")
+    print(
+        "Interpolation for a million coordinates should be"
+        " in the order of seconds. Anything significantly less "
+        "is considered slow"
+    )
     import time
+
     arr = np.linspace(-4, 4, 6000)
     arr = np.sin(arr)
-    arr = arr.reshape(10, 20, 30).astype('float32')
+    arr = arr.reshape(10, 20, 30).astype("float32")
     gridn = 1
     for c in range(8):
-        R1 = np.linspace(4., 5., gridn).astype('float32')
-        C1 = np.linspace(11., 12., gridn).astype('float32')
-        S1 = np.linspace(15., 16., gridn).astype('float32')
+        R1 = np.linspace(4.0, 5.0, gridn).astype("float32")
+        C1 = np.linspace(11.0, 12.0, gridn).astype("float32")
+        S1 = np.linspace(15.0, 16.0, gridn).astype("float32")
         tic = time.time()
         v1 = interp3(arr, R1, C1, S1)
         if gridn == 10 or gridn == 1:
             print(v1)
         toc = time.time()
-        print("1 followed by %d zeros" % c, "|", gridn, "|", \
-               toc - tic, "seconds")
+        print("1 followed by %d zeros" % c, "|", gridn, "|", toc - tic, "seconds")
         gridn = gridn * 10
diff --git a/setup.py b/setup.py
index 233492d..3e3a768 100644
--- a/setup.py
+++ b/setup.py
@@ -1,52 +1,34 @@
-from numpy.distutils.core import setup, Extension
-from numpy.distutils.misc_util import get_numpy_include_dirs 
+from setuptools import setup
 import os, sys
 
-mods = ['gradunwarp.core.coeffs', 'gradunwarp.core.globals', 
-        'gradunwarp.core.__init__', 'gradunwarp.__init__',
-        'gradunwarp.core.utils',
-        'gradunwarp.core.unwarp_resample',
-        'gradunwarp.core.gradient_unwarp',
-        'gradunwarp.core.tests.test_utils',
-       ]
+mods = [
+    "gradunwarp.core.coeffs",
+    "gradunwarp.core.globals",
+    "gradunwarp.core.__init__",
+    "gradunwarp.__init__",
+    "gradunwarp.core.utils",
+    "gradunwarp.core.unwarp_resample",
+    "gradunwarp.core.gradient_unwarp",
+]
 
-dats = [('gradunwarp/core/', ['gradunwarp/core/interp3_ext.c']),
-        ('gradunwarp/core/', ['gradunwarp/core/legendre_ext.c']),
-        ('gradunwarp/core/', ['gradunwarp/core/transform_coordinates_ext.c']),
-       ]
+scripts_cmd = [
+    "gradunwarp/core/gradient_unwarp.py",
+]
 
-# to build the C extension interp3_ext.c
-ext1 = Extension('gradunwarp.core.interp3_ext',
-                 include_dirs = get_numpy_include_dirs(),
-                 sources = ['gradunwarp/core/interp3_ext.c'],
-                 extra_compile_args=['-O3'])
-# to build the C extension legendre_ext.c
-ext2 = Extension('gradunwarp.core.legendre_ext',
-                 include_dirs = get_numpy_include_dirs(),
-                 sources = ['gradunwarp/core/legendre_ext.c'],
-                 extra_compile_args=['-O3'])
-# to build the C extension transform_coordinates_ext.c
-ext3 = Extension('gradunwarp.core.transform_coordinates_ext',
-                 include_dirs = get_numpy_include_dirs(),
-                 sources = ['gradunwarp/core/transform_coordinates_ext.c'],
-                 extra_compile_args=['-O3'])
 
-scripts_cmd = ['gradunwarp/core/gradient_unwarp.py',]
-        
-
-def configuration(parent_package='', top_path=None):
+def configuration(parent_package="", top_path=None):
     from numpy.distutils.misc_util import Configuration
-    config = Configuration('',parent_package,top_path)
-    config.add_data_files ( *dats )
+
+    config = Configuration("", parent_package, top_path)
     return config
 
-setup(name='gradunwarp',
-      version = '1.2.1+HCP',
-      description = 'HCP version of Gradient Unwarping Package for Python/Numpy',
-      author = 'Human Connectome Project',
-      py_modules  = mods,
-      ext_modules = [ext1, ext2, ext3],
-      scripts = scripts_cmd,
-      configuration=configuration,
-     )
 
+setup(
+    name="gradunwarp",
+    version="1.3.0+HCP",
+    description="HCP version of Gradient Unwarping Package for Python/Numpy",
+    author="Human Connectome Project",
+    py_modules=mods,
+    scripts=scripts_cmd,
+    configuration=configuration,
+)