diff --git a/gradunwarp/core/coeffs.py b/gradunwarp/core/coeffs.py index efa9a77..9ce6700 100644 --- a/gradunwarp/core/coeffs.py +++ b/gradunwarp/core/coeffs.py @@ -46,7 +46,7 @@ def coef_file_parse(cfile, txt_var_map): ''' # 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_re = re.compile(r'^[^\#]') # regex for first character not a '#' coef_file = open(cfile, 'r') for line in coef_file.readlines(): if coef_re.match(line): @@ -78,137 +78,80 @@ def get_siemens_coef(cfile): # allegra is slightly different 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)) - az = np.zeros((coef_array_sz, coef_array_sz)) - 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 = create_txt_var_map(coef_array_sz) coef_file_parse(cfile, txt_var_map) - return Coeffs(ax, ay, az, bx, by, bz, R0_m) + return Coeffs( + txt_var_map['alpha_y'], + txt_var_map['alpha_z'], + txt_var_map['alpha_x'], + txt_var_map['beta_y'], + txt_var_map['beta_x'], + txt_var_map['beta_z'], + R0_m) def get_ge_coef(cfile): ''' 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 = create_txt_var_map(coef_array_sz) coef_file_parse(cfile, txt_var_map) - return Coeffs(ax, ay, az, bx, by, bz, R0_m) + return Coeffs( + txt_var_map['alpha_y'], + txt_var_map['alpha_z'], + txt_var_map['alpha_x'], + txt_var_map['beta_y'], + txt_var_map['beta_x'], + txt_var_map['beta_z'], + R0_m) def grad_file_parse(gfile, txt_var_map): - ''' a separate function because GE and Siemens .coef files - have similar structure - - modifies txt_var_map in place - ''' - gf = open(gfile, 'r') - line = next(gf) - # skip the comments - while not line.startswith('#*] END:'): - line = next(gf) - - # get R0 - line = next(gf) - line = next(gf) - line = next(gf) - R0_m = float(line.strip().split()[0]) - - # go to the data - line = next(gf) - line = next(gf) - line = next(gf) - line = next(gf) - line = next(gf) - line = next(gf) - line = next(gf) - xmax = 0 ymax = 0 - - while 1: - 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) - 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]) - try: - line = next(gf) - except StopIteration: - break - - # just return R0_m but also txt_var_map is returned + with open(gfile, 'r') as gf: + for line in gf: + re_search = re.search(r"(?P\d+)[\s]+(?P[AB])[\s]*\(\s*(?P\d+),\s*(?P\d+)\)\s+(?P[\-]?\d+\.\d+)\s+(?P[xyz])", line) + if re_search: + re_res = re_search.groupdict() + alphabeta = 'alpha' if re_res['aorb'] == 'A' else 'beta' + x, y = int(re_res['x']),int(re_res['y']) + field = "%s_%s"%(alphabeta, re_res['axis']) + txt_var_map[field][x, y] = float(re_res['spectrum']) + xmax, ymax = max(x, xmax), max(y, ymax) + else: + re_search = re.search(r"(?P\d+\.\d+) m = R0", line) + if re_search: + R0_m = float(re_search.group('R0')) return R0_m, (xmax, ymax) +def create_txt_var_map(coef_array_sz): + txt_var_map = {} + for ab in ['alpha', 'beta']: + for axis in 'xyz': + txt_var_map['%s_%s'%(ab,axis)] = np.zeros((coef_array_sz,)*2) + return txt_var_map + def get_siemens_grad(gfile): ''' Parse the siemens .grad file ''' - coef_array_sz = siemens_cas # allegra is slightly different - 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)) - az = np.zeros((coef_array_sz, coef_array_sz)) - 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} + coef_array_sz = 15 if gfile.startswith('coef_AC44') else siemens_cas - 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] - - return Coeffs(ax, ay, az, bx, by, bz, R0_m) + txt_var_map = create_txt_var_map(coef_array_sz) + print(txt_var_map) + R0_m, max_ind = grad_file_parse(gfile, txt_var_map) + ind = max(max_ind) + 1 + + return Coeffs( + txt_var_map['alpha_x'][:ind, :ind], + txt_var_map['alpha_y'][:ind, :ind], + txt_var_map['alpha_z'][:ind, :ind], + txt_var_map['beta_x'][:ind, :ind], + txt_var_map['beta_y'][:ind, :ind], + txt_var_map['beta_z'], + R0_m) 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 - -# 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 - -/* 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= 3) - { - k -= 2; - f *= k; - } - return f; -} diff --git a/gradunwarp/core/tests/data/gradunwarp_coeffs.grad b/gradunwarp/core/tests/data/gradunwarp_coeffs.grad new file mode 100644 index 0000000..2d4da85 --- /dev/null +++ b/gradunwarp/core/tests/data/gradunwarp_coeffs.grad @@ -0,0 +1,37 @@ +#*[ Script ****************************************************************\ +# +# Name : dummy.grad +# +# Date : 0970-01-01 +# +# Author : Mario Bros +# +# Language : +# +# Description : Defines Legendre coefficients in spherical harmonics for +# dummy Gradient Coil +# +#****************************************************************************/ + +#*] END: */ +dummy, dummy , Gx,y,z = 80/80 mT/m +win_low = 0, win_high = 0, win_algo = 0, win_dummy = 2; +0.25 m = R0, lnorm = 4? A(1,0) = B(1,1) = A(1,1) = 0; +0 = CoSyMode, + + + +NO. TYPE SPECTRUM AXIS + + 1 A( 3, 0) -0.0100000 z + 2 A( 5, 0) 0.02000000 z + 3 A( 7, 0) -0.0300000 z + 4 A( 9, 0) 0.0400000 z +101 A( 3, 1) 0.0100000 x +102 A( 3, 3) -0.02000000 x +103 A( 5, 1) 0.03000000 x +104 A( 5, 3) -0.04000000 x +201 B( 3, 1) 0.01000000 y +202 B( 3, 3) -0.02000000 y +203 B( 5, 1) 0.03000000 y +204 B( 5, 3) -0.04000000 y diff --git a/gradunwarp/core/tests/data/siemens_B_output.npz b/gradunwarp/core/tests/data/siemens_B_output.npz new file mode 100644 index 0000000..558bde7 Binary files /dev/null and b/gradunwarp/core/tests/data/siemens_B_output.npz differ diff --git a/gradunwarp/core/tests/test_regression.py b/gradunwarp/core/tests/test_regression.py new file mode 100644 index 0000000..9854a4c --- /dev/null +++ b/gradunwarp/core/tests/test_regression.py @@ -0,0 +1,28 @@ +import os +DATA_DIR = os.path.join(os.path.dirname(__file__), 'data') + +import numpy as np +from numpy.testing import assert_array_almost_equal, assert_allclose + +from gradunwarp.core import coeffs, utils +from gradunwarp.core.unwarp_resample import siemens_B, cart2sph + +def test_siemens_B(): + gradfile = os.path.join(DATA_DIR, 'gradunwarp_coeffs.grad') + + siemens_coeffs = coeffs.get_coefficients('siemens', gradfile) + R0 = siemens_coeffs.R0_m * 1000 + + vec = np.linspace(-300, 300, 60, dtype=np.float32) + x, y ,z = np.meshgrid(vec, vec, vec) + r, cos_theta, theta, phi = cart2sph(x, y, z) + + ref_b = np.load(os.path.join(DATA_DIR, 'siemens_B_output.npz')) + + for d in 'xyz': + alpha_d = getattr(siemens_coeffs, "alpha_%s"%d) + beta_d = getattr(siemens_coeffs, "beta_%s"%d) + bd = siemens_B(alpha_d, beta_d, r, cos_theta, theta, phi, R0) + + # changes in legendre function is causing differences at 5th decimal + assert_array_almost_equal(ref_b["b%s"%d], bd, decimal=5) 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 - -/* 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 87b01bc..4782800 100644 --- a/gradunwarp/core/unwarp_resample.py +++ b/gradunwarp/core/unwarp_resample.py @@ -19,6 +19,7 @@ from .globals import siemens_max_det import nibabel as nib import subprocess +import scipy.special #np.seterr(all='raise') @@ -49,7 +50,7 @@ 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. ''' @@ -72,10 +73,10 @@ def eval_spharm_grid(self, vendor, coeffs): fovmax = fovmax * 1000. # the grid in meters. this is needed for spherical harmonics vec = np.linspace(fovmin, fovmax, numpoints) - gvx, gvy, gvz = utils.meshgrid(vec, vec, vec) + gvx, gvy, gvz = np.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], @@ -85,16 +86,12 @@ def eval_spharm_grid(self, vendor, coeffs): # 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') gvxyz = CV(gvx, gvy, gvz) _dv, _dxyz = eval_spherical_harmonics(coeffs, vendor, gvxyz) - + return CV(_dv.x, _dv.y, _dv.z), g_xyz2rcs @@ -107,7 +104,7 @@ def run(self): if self.warp: self.polarity = -1. - # 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 @@ -118,14 +115,6 @@ def run(self): 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 @@ -146,7 +135,7 @@ def run(self): # through rotation and scaling, where any mirror # will impose a negation ones = CV(1., 1., 1.) - dxyz = utils.transform_coordinates_old(ones, r_rcs2lai) + dxyz = utils.transform_coordinates(ones, r_rcs2lai) # do the nonlinear unwarp if self.vendor == 'siemens': @@ -178,25 +167,25 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh 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) - vjacdet_lpsw = np.zeros((nr, nc), dtype=np.float32) + jim2 = np.empty((nr, nc), dtype=np.float32) + vjacdet_lpsw = np.empty((nr, nc), dtype=np.float32) if dxyz == 0: vjacdet_lps = 1 else: vjacdet_lps = eval_siemens_jacobian_mult(dv, dxyz) - # essentially pre-allocating everything - out = np.zeros((nr, nc, ns), dtype=np.float32) - fullWarp = np.zeros((nr, nc, ns, 3), dtype=np.float32) + # essentially pre-allocating everything + out = np.empty((nr, nc, ns), dtype=np.float32) + fullWarp = np.empty((nr, nc, ns, 3), dtype=np.float32) - vjacout = np.zeros((nr, nc, ns), dtype=np.float32) - im2 = np.zeros((nr, nc), dtype=np.float32) - dvx = np.zeros((nr, nc), dtype=np.float32) - dvy = np.zeros((nr, nc), dtype=np.float32) - dvz = np.zeros((nr, nc), dtype=np.float32) - im_ = np.zeros((nr, nc), dtype=np.float32) + vjacout = np.empty((nr, nc, ns), dtype=np.float32) + im2 = np.empty((nr, nc), dtype=np.float32) + dvx = np.empty((nr, nc), dtype=np.float32) + dvy = np.empty((nr, nc), dtype=np.float32) + dvz = np.empty((nr, nc), dtype=np.float32) + im_ = np.empty((nr, nc), dtype=np.float32) # init jacobian temp image - vc, vr = utils.meshgrid(np.arange(nc), np.arange(nr)) + vc, vr = np.meshgrid(np.arange(nc), np.arange(nr)) # Compute transform to map the internal voxel coordinates to FSL scaled mm coordinates img = nib.load(self.name) @@ -215,7 +204,7 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh [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): @@ -225,37 +214,25 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh print(s+1, end=' ') else: print('.', end=' ') - - # hopefully, free memory - gc.collect() - # init to 0 - dvx.fill(0.) - dvy.fill(0.) - dvz.fill(0.) - im_.fill(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.map_coordinates(dv.x, + vrcsg, + output=dvx, + order=self.order) + ndimage.map_coordinates(dv.y, + vrcsg, + output=dvy, + order=self.order) + ndimage.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) @@ -268,36 +245,28 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh #im_ = utils.interp3(self.vol, vrcsw.x, vrcsw.y, vrcsw.z) - ndimage.interpolation.map_coordinates(self.vol, - vrcsw, - output=im_, - order=self.order) + ndimage.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. + np.nan_to_num(im_, copy=False) + im_[np.isinf(im_)] = 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") - # Multiply the intensity with the Jacobian det, if needed if not self.nojac: - vjacdet_lpsw.fill(0.) jim2.fill(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. + ndimage.map_coordinates(vjacdet_lps, + vrcsg, + output=vjacdet_lpsw, + order=self.order) + np.nan_to_num(vjacdet_lpsw, copy=False) + vjacdet_lpsw[np.isinf(vjacdet_lpsw)] = 0 jim2[vr, vc] = vjacdet_lpsw im2 = im2 * jim2 vjacout[..., s] = jim2 @@ -308,8 +277,8 @@ def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai_noh out[..., s] = im2 print() - - img=nib.Nifti1Image(fullWarp,self.m_rcs2ras) + + img=nib.Nifti1Image(fullWarp, self.m_rcs2ras) nib.save(img,"fullWarp_abs.nii.gz") # return image and the jacobian return out, vjacout @@ -371,72 +340,66 @@ def eval_spherical_harmonics(coeffs, vendor, vxyz): x, y, z = vxyz - #pdb.set_trace() - # log.info('calculating displacements (mm) ' - # 'using spherical harmonics coeffcients...') + r, cos_theta, theta, phi = cart2sph(x, y, z) + if vendor == 'siemens': log.info('along x...') - bx = siemens_B(coeffs.alpha_x, coeffs.beta_x, x, y, z, R0) + bx = siemens_B(coeffs.alpha_x, coeffs.beta_x, r, cos_theta, theta, phi, R0) log.info('along y...') - by = siemens_B(coeffs.alpha_y, coeffs.beta_y, x, y, z, R0) + by = siemens_B(coeffs.alpha_y, coeffs.beta_y, r, cos_theta, theta, phi, R0) log.info('along z...') - bz = siemens_B(coeffs.alpha_z, coeffs.beta_z, x, y, z, R0) + bz = siemens_B(coeffs.alpha_z, coeffs.beta_z, r, cos_theta, theta, phi, R0) else: # GE log.info('along x...') - bx = ge_D(coeffs.alpha_x, coeffs.beta_x, x, y, z) + bx = ge_D(coeffs.alpha_x, coeffs.beta_x, r, cos_theta, theta, phi) log.info('along y...') - by = ge_D(coeffs.alpha_y, coeffs.beta_y, x, y, z) + by = ge_D(coeffs.alpha_y, coeffs.beta_y, r, cos_theta, theta, phi) 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) + bz = ge_D(coeffs.alpha_z, coeffs.beta_z, r, cos_theta, theta, phi) return CV(bx * R0, by * R0, bz * R0), CV(x, y, z) -#@profile -def siemens_B(alpha, beta, x1, y1, z1, R0): - ''' Calculate displacement field from Siemens coefficients - ''' - nmax = alpha.shape[0] - 1 +def cart2sph(x1, y1, z1): x1 = x1 + 0.0001 # hack to avoid singularities at R=0 - # convert to spherical coordinates r = np.sqrt(x1 * x1 + y1 * y1 + z1 * z1) - theta = np.arccos(z1 / r) + cosine_theta = z1 / r + theta = np.arccos(cosine_theta) phi = np.arctan2(y1 / r, x1 / r) - b = np.zeros(x1.shape) + return r, cosine_theta, theta, phi + +def siemens_B(alpha, beta, r, cosine_theta, theta, phi, R0): + ''' Calculate displacement field from Siemens coefficients + ''' + nmax = alpha.shape[0] - 1 + + b = np.zeros(theta.shape) for n in range(0, nmax + 1): 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)) - normfact = 1 + _p = scipy.special.lpmv(m, n, cosine_theta) # 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 + _p *= normfact + b += f * _p * f2 return b - -def ge_D(alpha, beta, x1, y1, z1): +# For consistency with GE papers, use theta & phi -> phi & theta +def ge_D(alpha, beta, r, cosine_phi, phi, theta, z1): ''' 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 - r = np.sqrt(x1 * x1 + y1 * y1 + z1 * z1) - # For consistency with GE papers, use theta & phi -> phi & theta - phi = np.arccos(z1 / r) - theta = np.arctan2(y1 / r, x1 / r) r = r * 100.0 # GE wants cm, so meters -> cm - d = np.zeros(x1.shape) + d = np.zeros(theta.shape) for n in range(0, nmax + 1): # So GE uses the usual unnormalized legendre polys. @@ -444,7 +407,7 @@ def ge_D(alpha, beta, x1, y1, z1): 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 + _p = scipy.special.lpmv(m, n, cosine_phi) + d += f * _p * f2 d = d / 100.0 # cm back to meters return d diff --git a/gradunwarp/core/utils.py b/gradunwarp/core/utils.py index 42bde47..7379640 100644 --- a/gradunwarp/core/utils.py +++ b/gradunwarp/core/utils.py @@ -18,42 +18,11 @@ # cv = CoordsVector(x=x, y=y, z=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 -# are costly for huge matrices -def transform_coordinates_old(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 - B1 = A1 * M[0, 0] + A2 * M[0, 1] + A3 * M[0, 2] + M[0, 3] - B2 = A1 * M[1, 0] + A2 * M[1, 1] + A3 * M[1, 2] + M[1, 3] - B3 = A1 * M[2, 0] + A2 * M[2, 1] + A3 * M[2, 2] + M[2, 3] - return CoordsVector(B1, B2, B3) +from nibabel.affines import apply_affine def transform_coordinates(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 - A1 = A1.astype(np.float32) - A2 = A2.astype(np.float32) - A3 = A3.astype(np.float32) - M = M.astype(np.float32) - 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.') - - B1, B2, B3 = _transform_coordinates(A1, A2, A3, M) - return CoordsVector(B1, B2, B3) - + transformed = apply_affine(M, np.stack(A).T).T + return CoordsVector(*transformed) def get_vol_affine(infile): try: @@ -64,264 +33,4 @@ def get_vol_affine(infile): nibimage = nib.load(infile) return np.asanyarray(nibimage.dataobj), nibimage.affine - -# memoized factorial -class Memoize: - def __init__(self, f): - self.f = f - self.memo = {} - - def __call__(self, *args): - if not args in self.memo: - self.memo[args] = self.f(*args) - return self.memo[args] - -#factorial = Memoize(math.factorial) factorial = math.factorial - - -### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -# -# This is taken from a numpy ticket which adds N-d support to meshgrid -# URL : http://projects.scipy.org/numpy/ticket/966 -# License : http://docs.scipy.org/doc/numpy/license.html -# -### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ### ## -def meshgrid(*xi, **kwargs): - """ - Return coordinate matrices from one or more coordinate vectors. - - Make N-D coordinate arrays for vectorized evaluations of - N-D scalar/vector fields over N-D grids, given - one-dimensional coordinate arrays x1, x2,..., xn. - - Parameters - ---------- - x1, x2,..., xn : array_like - 1-D arrays representing the coordinates of a grid. - indexing : 'xy' or 'ij' (optional) - cartesian ('xy', default) or matrix ('ij') indexing of output - sparse : True or False (default) (optional) - If True a sparse grid is returned in order to conserve memory. - copy : True (default) or False (optional) - If False a view into the original arrays are returned in order to - conserve memory - - Returns - ------- - X1, X2,..., XN : ndarray - For vectors `x1`, `x2`,..., 'xn' with lengths ``Ni=len(xi)`` , - return ``(N1, N2, N3,...Nn)`` shaped arrays if indexing='ij' - or ``(N2, N1, N3,...Nn)`` shaped arrays if indexing='xy' - with the elements of `xi` repeated to fill the matrix along - the first dimension for `x1`, the second for `x2` and so on. - - See Also - -------- - index_tricks.mgrid : Construct a multi-dimensional "meshgrid" - using indexing notation. - index_tricks.ogrid : Construct an open multi-dimensional "meshgrid" - using indexing notation. - - Examples - -------- - >>> x = np.linspace(0,1,3) # coordinates along x axis - >>> y = np.linspace(0,1,2) # coordinates along y axis - >>> xv, yv = meshgrid(x,y) # extend x and y for a 2D xy grid - >>> xv - array([[ 0. , 0.5, 1. ], - [ 0. , 0.5, 1. ]]) - >>> yv - array([[ 0., 0., 0.], - [ 1., 1., 1.]]) - >>> xv, yv = meshgrid(x,y, sparse=True) # make sparse output arrays - >>> xv - array([[ 0. , 0.5, 1. ]]) - >>> yv - array([[ 0.], - [ 1.]]) - - >>> meshgrid(x,y,sparse=True,indexing='ij') # change to matrix indexing - [array([[ 0. ], - [ 0.5], - [ 1. ]]), array([[ 0., 1.]])] - >>> meshgrid(x,y,indexing='ij') - [array([[ 0. , 0. ], - [ 0.5, 0.5], - [ 1. , 1. ]]), array([[ 0., 1.], - [ 0., 1.], - [ 0., 1.]])] - - >>> meshgrid(0,1,5) # just a 3D point - [array([[[0]]]), array([[[1]]]), array([[[5]]])] - >>> map(np.squeeze,meshgrid(0,1,5)) # just a 3D point - [array(0), array(1), array(5)] - >>> meshgrid(3) - array([3]) - >>> meshgrid(y) # 1D grid; y is just returned - array([ 0., 1.]) - - `meshgrid` is very useful to evaluate functions on a grid. - - >>> x = np.arange(-5, 5, 0.1) - >>> y = np.arange(-5, 5, 0.1) - >>> xx, yy = meshgrid(x, y, sparse=True) - >>> z = np.sin(xx**2+yy**2)/(xx**2+yy**2) - """ - 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)') - - 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)] - - shape = [x.size for x in output] - - if indexing == 'xy': - # switch first and second axis - 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: - if copy: - return [x.copy() for x in output] - else: - return output - else: - # Return the full N-D matrix (not only the 1-D vector) - if copy: - mult_fact = np.ones(shape, dtype=int) - return [x * mult_fact for x in output] - else: - return np.broadcast_arrays(*output) - - -def ndgrid(*args, **kwargs): - """ - Same as calling meshgrid with indexing='ij' (see meshgrid for - documentation). - """ - kwargs['indexing'] = 'ij' - return meshgrid(*args, **kwargs) - - -def odd_factorial(k): - f = k - while k >= 3: - k -= 2 - f *= k - return f - - -# From the scipy ticket -# http://projects.scipy.org/scipy/attachment/ticket/1296/assoc_legendre.py -def legendre_old(nu, mu, x): - """Compute the associated Legendre polynomial with degree nu and order mu. - This function uses the recursion formula in the degree nu. - (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)) - # if abs(x) > 1: - # raise ValueError('require -1 <= x <= 1, but x=%f', x) - - # Compute the initial term in the recursion. - if mu == 0: - p_nu = 1.0 - else: - s = 1 - if mu & 1: - s = -1 - z = sqrt(1 - x ** 2) - p_nu = s * odd_factorial(2 * mu - 1) * z ** mu - - if mu == nu: - return p_nu - - # Compute the next term in the recursion. - p_nu_prev = p_nu - p_nu = x * (2 * mu + 1) * p_nu - - if nu == mu + 1: - return p_nu - - # Iterate the recursion relation. - for n in range(mu + 2, nu + 1): - result = (x * (2 * n - 1) * p_nu - (n + mu - 1) * p_nu_prev) / (n - mu) - p_nu_prev = p_nu - p_nu = result - - return result - - -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.') - 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.') - vol = vol.astype(np.float32) - vol = np.ascontiguousarray(vol) - R = R.astype(np.float32) - C = C.astype(np.float32) - S = S.astype(np.float32) - assert vol.ndim == 3 - assert C.ndim == R.ndim == S.ndim - assert C.shape == R.shape == S.shape - V = _interp3(vol, R, C, S) - return V - - -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') - import time - arr = np.linspace(-4, 4, 6000) - arr = np.sin(arr) - 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') - 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") - gridn = gridn * 10 diff --git a/setup.cfg b/setup.cfg index cfa9efa..39fff7d 100644 --- a/setup.cfg +++ b/setup.cfg @@ -28,3 +28,6 @@ install_requires = numpy>=1.26.0b1; python_version > '3.11' numpy; python_version <= '3.11' scipy + +[options.package_data] +gradunwarp.core.tests = data/* diff --git a/setup.py b/setup.py index 21ffd0e..9e554bb 100644 --- a/setup.py +++ b/setup.py @@ -3,19 +3,8 @@ Most configuration is now in pyproject.toml. This file configures extensions and a legacy script. """ -from setuptools import setup, Extension -from numpy import get_include - +from setuptools import setup setup( - ext_modules=[ - Extension( - 'gradunwarp.core.{}_ext'.format(mod), - include_dirs=[get_include()], - sources=['gradunwarp/core/{}_ext.c'.format(mod)], - extra_compile_args=['-O3'], - ) - for mod in ('interp3', 'legendre', 'transform_coordinates') - ], scripts=['gradunwarp/core/gradient_unwarp.py'], )