From a49d564ca9edbacef17057489ee2610614d66196 Mon Sep 17 00:00:00 2001 From: Daniel Rose Date: Mon, 12 Oct 2020 11:08:21 +0200 Subject: [PATCH 1/5] added and moved files around to make a python package out of this repo --- .gitignore | 140 ++++++++++++++++++ example_scripts/run_SM_fit_example.py | 3 +- example_scripts/validate_formulas_example.py | 4 +- linear_planar/__init__.py | 0 b02b0_2mm.cnf => linear_planar/b02b0_2mm.cnf | 0 .../calc_jacobian.py | 0 debias.py => linear_planar/debias.py | 0 .../eddy_nlgc_mp.py | 0 gs.py => linear_planar/gs.py | 0 .../mean_std_4th.py | 0 .../smt_lin_pla.py | 0 .../split_shells_and_types.py | 0 setup.py | 27 ++++ 13 files changed, 170 insertions(+), 4 deletions(-) create mode 100755 .gitignore create mode 100644 linear_planar/__init__.py rename b02b0_2mm.cnf => linear_planar/b02b0_2mm.cnf (100%) rename calc_jacobian.py => linear_planar/calc_jacobian.py (100%) rename debias.py => linear_planar/debias.py (100%) rename eddy_nlgc_mp.py => linear_planar/eddy_nlgc_mp.py (100%) rename gs.py => linear_planar/gs.py (100%) rename mean_std_4th.py => linear_planar/mean_std_4th.py (100%) rename smt_lin_pla.py => linear_planar/smt_lin_pla.py (100%) rename split_shells_and_types.py => linear_planar/split_shells_and_types.py (100%) create mode 100755 setup.py diff --git a/.gitignore b/.gitignore new file mode 100755 index 0000000..7d3a2df --- /dev/null +++ b/.gitignore @@ -0,0 +1,140 @@ +# Created by .ignore support plugin (hsz.mobi) +### Example user template template +### Example user template + +# IntelliJ project files +.idea +*.iml +out +gen +### Python template +# Byte-compiled / optimized / DLL files +__pycache__/ +*.py[cod] +*$py.class + +# C extensions +*.so + +# Distribution / packaging +.Python +build/ +develop-eggs/ +dist/ +downloads/ +eggs/ +.eggs/ +lib/ +lib64/ +parts/ +sdist/ +var/ +wheels/ +pip-wheel-metadata/ +share/python-wheels/ +*.egg-info/ +.installed.cfg +*.egg +MANIFEST + +# PyInstaller +# Usually these files are written by a python script from a template +# before PyInstaller builds the exe, so as to inject date/other infos into it. +*.manifest +*.spec + +# Installer logs +pip-log.txt +pip-delete-this-directory.txt + +# Unit test / coverage reports +htmlcov/ +.tox/ +.nox/ +.coverage +.coverage.* +.cache +nosetests.xml +coverage.xml +*.cover +*.py,cover +.hypothesis/ +.pytest_cache/ + +# Translations +*.mo +*.pot + +# Django stuff: +*.log +local_settings.py +db.sqlite3 +db.sqlite3-journal + +# Flask stuff: +instance/ +.webassets-cache + +# Scrapy stuff: +.scrapy + +# Sphinx documentation +docs/_build/ + +# PyBuilder +target/ + +# Jupyter Notebook +.ipynb_checkpoints + +# IPython +profile_default/ +ipython_config.py + +# pyenv +.python-version + +# pipenv +# According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control. +# However, in case of collaboration, if having platform-specific dependencies or dependencies +# having no cross-platform support, pipenv may install dependencies that don't work, or not +# install all needed dependencies. +#Pipfile.lock + +# PEP 582; used by e.g. github.com/David-OConnor/pyflow +__pypackages__/ + +# Celery stuff +celerybeat-schedule +celerybeat.pid + +# SageMath parsed files +*.sage.py + +# Environments +.env +.venv +env/ +venv/ +ENV/ +env.bak/ +venv.bak/ + +# Spyder project settings +.spyderproject +.spyproject + +# Rope project settings +.ropeproject + +# mkdocs documentation +/site + +# mypy +.mypy_cache/ +.dmypy.json +dmypy.json + +# Pyre type checker +.pyre/ + diff --git a/example_scripts/run_SM_fit_example.py b/example_scripts/run_SM_fit_example.py index 6dbbaf0..30f2cdb 100644 --- a/example_scripts/run_SM_fit_example.py +++ b/example_scripts/run_SM_fit_example.py @@ -3,11 +3,10 @@ import scipy.optimize as opt # import pylab as pl -from smt_lin_pla import * +from linear_planar.smt_lin_pla import * from multiprocess import Pool -import sys import os from time import time diff --git a/example_scripts/validate_formulas_example.py b/example_scripts/validate_formulas_example.py index 706fe06..f0eab84 100644 --- a/example_scripts/validate_formulas_example.py +++ b/example_scripts/validate_formulas_example.py @@ -1,8 +1,8 @@ import numpy as np -from gs import findLinIndepRandomRot, gramSchmidt3 +from linear_planar.gs import findLinIndepRandomRot, gramSchmidt3 from dipy.core.sphere import HemiSphere, disperse_charges # from sphHist import plotScatter3 -from smt_lin_pla import * +from linear_planar.smt_lin_pla import * ## we want to generate a set of well distribution point on the sphere diff --git a/linear_planar/__init__.py b/linear_planar/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/b02b0_2mm.cnf b/linear_planar/b02b0_2mm.cnf similarity index 100% rename from b02b0_2mm.cnf rename to linear_planar/b02b0_2mm.cnf diff --git a/calc_jacobian.py b/linear_planar/calc_jacobian.py similarity index 100% rename from calc_jacobian.py rename to linear_planar/calc_jacobian.py diff --git a/debias.py b/linear_planar/debias.py similarity index 100% rename from debias.py rename to linear_planar/debias.py diff --git a/eddy_nlgc_mp.py b/linear_planar/eddy_nlgc_mp.py similarity index 100% rename from eddy_nlgc_mp.py rename to linear_planar/eddy_nlgc_mp.py diff --git a/gs.py b/linear_planar/gs.py similarity index 100% rename from gs.py rename to linear_planar/gs.py diff --git a/mean_std_4th.py b/linear_planar/mean_std_4th.py similarity index 100% rename from mean_std_4th.py rename to linear_planar/mean_std_4th.py diff --git a/smt_lin_pla.py b/linear_planar/smt_lin_pla.py similarity index 100% rename from smt_lin_pla.py rename to linear_planar/smt_lin_pla.py diff --git a/split_shells_and_types.py b/linear_planar/split_shells_and_types.py similarity index 100% rename from split_shells_and_types.py rename to linear_planar/split_shells_and_types.py diff --git a/setup.py b/setup.py new file mode 100755 index 0000000..d85f8f7 --- /dev/null +++ b/setup.py @@ -0,0 +1,27 @@ +#! /usr/bin/env python +# -*- coding: utf-8 -*- + +"""Python install script. + +Install with `python setup.py install` or `python setup.py develop` for developer mode.""" + +from setuptools import setup, find_packages + +REQUIREMENTS = ["numpy", + "nibabel", + "scipy"] + +setup( + name="linear_planar", + description="Preprocessing and fitting scripts for double PGSE data", + version="0.1", + packages=find_packages(), + package_data={ + "": ["*.txt", "*.cnf", "*.sh"], + }, + install_requires=REQUIREMENTS, + python_requires=">=3.5", + url="https://github.com/blochchainlab/linear_planar" +) + + From 86677512acc682dbdff598ba009d9dacb28fcf9c Mon Sep 17 00:00:00 2001 From: Daniel Rose Date: Mon, 12 Oct 2020 11:11:06 +0200 Subject: [PATCH 2/5] fixed indents to 4 spaces instead of tabs and fixed some PEP8 compliance (auto-refomatting by PyCharm) --- example_scripts/run_SM_fit_example.py | 167 +++++++------------ example_scripts/validate_formulas_example.py | 101 ++++++----- linear_planar/calc_jacobian.py | 25 ++- linear_planar/eddy_nlgc_mp.py | 53 +++--- linear_planar/gs.py | 93 ++++++----- linear_planar/mean_std_4th.py | 20 +-- linear_planar/smt_lin_pla.py | 80 ++++----- linear_planar/split_shells_and_types.py | 107 ++++++------ 8 files changed, 297 insertions(+), 349 deletions(-) diff --git a/example_scripts/run_SM_fit_example.py b/example_scripts/run_SM_fit_example.py index 30f2cdb..9cea9d9 100644 --- a/example_scripts/run_SM_fit_example.py +++ b/example_scripts/run_SM_fit_example.py @@ -11,17 +11,9 @@ from time import time - - ## for multiprocess N_pool = 6 - - - - - - ## This part would be kind-of how to create b and bd from G1G2G3.txt instead of hardcoding # gtable = np.loadtxt('/home/raid2/paquette/work/genDirsPlanar/scheme/test_45_45_x3.txt') # ten = table2ten(gtable, bval=1.5) @@ -31,12 +23,6 @@ # scheme, scheme_norm, dir_type = bd_getDirs(ten) - - - - - - ## load data masterpath = '/some/path/' @@ -45,64 +31,56 @@ # this structure assumes that each bvalue and b-shape are save separately in already spherical averaged form data_path = [masterpath + 'PREPRO/v1_eddy_b1000_lin_SM.nii.gz', - masterpath + 'PREPRO/v1_eddy_b1000_pla_SM.nii.gz', - masterpath + 'PREPRO/v1_eddy_b2000_lin_SM.nii.gz', - masterpath + 'PREPRO/v1_eddy_b2000_pla_SM.nii.gz', - masterpath + 'PREPRO/v1_eddy_b4000_lin_SM.nii.gz', - masterpath + 'PREPRO/v1_eddy_b4000_pla_SM.nii.gz'] + masterpath + 'PREPRO/v1_eddy_b1000_pla_SM.nii.gz', + masterpath + 'PREPRO/v1_eddy_b2000_lin_SM.nii.gz', + masterpath + 'PREPRO/v1_eddy_b2000_pla_SM.nii.gz', + masterpath + 'PREPRO/v1_eddy_b4000_lin_SM.nii.gz', + masterpath + 'PREPRO/v1_eddy_b4000_pla_SM.nii.gz'] mask_path = masterpath + 'PREPRO/bet_init_b50_dil.nii.gz' savepath = masterpath + 'FIT/fit1/' - ## hardcoded b-tensor size and shape corresponding to the load order of data_path (in um^2/ms) b = np.array([1., 1., 2., 2., 4., 4.]) bd = np.array([1., -.5, 1., -.5, 1., -.5]) - if not os.path.isdir(savepath): - print('saving directory doesn\'t exist') - # return None + print('saving directory doesn\'t exist') +# return None b0_ima = nib.load(b0_path) b0_data = b0_ima.get_data() -data = np.concatenate([nib.load(path).get_data()[...,None] for path in data_path], axis=3) +data = np.concatenate([nib.load(path).get_data()[..., None] for path in data_path], axis=3) mask = nib.load(mask_path).get_data().astype(bool) print('data loaded') - - - ## data b0 normalization, clip outlier to [0, 1] -data_norm = np.clip(data, 0, np.inf) / np.clip(b0_data, 0, np.inf)[...,None] +data_norm = np.clip(data, 0, np.inf) / np.clip(b0_data, 0, np.inf)[..., None] data_norm[np.isnan(data_norm)] = 0 # clip outlier data_SM = np.clip(data_norm, 0, 1) print('data normalized') - - - ## defining objective function from signal function of smt_lin_pla (model2) ## note: if using model 1, f_csf is no longuer a parameter # single-voxel per-shell error function def func_vox(p, b, bd, datavox): - # unpack param - Dpar_in, Dpar_ex, Dperp_ex, f_in, f_ex, f_csf = p - # return model(param)-data for model2 - return sm_signal_model_2(b, bd, Dpar_in, Dpar_ex, Dperp_ex, f_in, f_ex, f_csf) - datavox + # unpack param + Dpar_in, Dpar_ex, Dperp_ex, f_in, f_ex, f_csf = p + # return model(param)-data for model2 + return sm_signal_model_2(b, bd, Dpar_in, Dpar_ex, Dperp_ex, f_in, f_ex, f_csf) - datavox -# single-voxel squared_error -def func_vox_ls(p, b, bd, datavox): - # compute single-voxel per-shell error function - f = func_vox(p, b, bd, datavox) - # return sum_{shell} (model(param)-data)^2 - return np.sum(f**2, axis=0) +# single-voxel squared_error +def func_vox_ls(p, b, bd, datavox): + # compute single-voxel per-shell error function + f = func_vox(p, b, bd, datavox) + # return sum_{shell} (model(param)-data)^2 + return np.sum(f ** 2, axis=0) ################################################# @@ -150,24 +128,22 @@ def func_vox_ls(p, b, bd, datavox): # 1 <= f_in + f_ex + f_csf <= 1 A = np.array([[1, 0, 0, 0, 0, 0], - [0, 1, 0, 0, 0, 0], - [0, 0, 1, 0, 0, 0], - [0, 0, 0, 1, 0, 0], - [0, 0, 0, 0, 1, 0], - [0, 0, 0, 0, 0, 1], - [0, 0, 0, 1, 1, 1]]) - -lb = np.array([0.1,0.1,0.1,0,0,0,1]) -ub = np.array([3,3,1.5,1,1,1,1]) + [0, 1, 0, 0, 0, 0], + [0, 0, 1, 0, 0, 0], + [0, 0, 0, 1, 0, 0], + [0, 0, 0, 0, 1, 0], + [0, 0, 0, 0, 0, 1], + [0, 0, 0, 1, 1, 1]]) + +lb = np.array([0.1, 0.1, 0.1, 0, 0, 0, 1]) +ub = np.array([3, 3, 1.5, 1, 1, 1, 1]) ################################################# - ## linear constraint object cons = opt.LinearConstraint(A, lb, ub) - ################################################# ## non-linear ratio constraint set @@ -179,20 +155,22 @@ def func_vox_ls(p, b, bd, datavox): def ratio(x): - # Dpar / Dperp - return x[1]/x[2] + # Dpar / Dperp + return x[1] / x[2] + # derivative of the constraint for all parameters def ratio_jac(x): - # d/dDin = 0 - # d/dDpar = 1/Dperp - d1 = 1 / x[2] - # d/dPerp = -Dpar / Dperp^2 - d2 = -x[1] / x[2]**2 - # d/dfin = 0 - # d/dfex = 0 - # d/dfw = 0 - return np.array([0,d1,d2,0,0,0]) + # d/dDin = 0 + # d/dDpar = 1/Dperp + d1 = 1 / x[2] + # d/dPerp = -Dpar / Dperp^2 + d2 = -x[1] / x[2] ** 2 + # d/dfin = 0 + # d/dfex = 0 + # d/dfw = 0 + return np.array([0, d1, d2, 0, 0, 0]) + ## some FA to ratio conversion, computed elsewhere # # FA ~ 0.2425 @@ -210,31 +188,24 @@ def ratio_jac(x): fa_cons = opt.NonlinearConstraint(ratio, np.array([ratio_min]), np.array([ratio_max]), jac=ratio_jac, hess=opt.BFGS()) - - - - ## uFA for this type of models def microFA(d1, d2, d3, fin, fex, fcsf): - # if (fin < 0) or (fex < 0) or (fin + fex > 1): - # return np.nan + # if (fin < 0) or (fex < 0) or (fin + fex > 1): + # return np.nan - microAx = fin*d1 + fex*d2 + fcsf*3 - microRad = fex*d3 + fcsf*3 - microMean = (microAx + 2*microRad) / 3. - microFA = np.sqrt((3/2.)*((microAx - microMean)**2 + 2*(microRad - microMean)**2) / (microAx**2 + 2*microRad**2)) - return microFA + microAx = fin * d1 + fex * d2 + fcsf * 3 + microRad = fex * d3 + fcsf * 3 + microMean = (microAx + 2 * microRad) / 3. + microFA = np.sqrt( + (3 / 2.) * ((microAx - microMean) ** 2 + 2 * (microRad - microMean) ** 2) / (microAx ** 2 + 2 * microRad ** 2)) + return microFA ## compartement FA def FA(Dpar, Dperp): - meanD = (Dpar+2*Dperp)/3. - return np.sqrt(3/2.) * np.sqrt((Dpar-meanD)**2 + 2*(Dperp-meanD)**2) / np.sqrt(Dpar**2 + 2*Dperp**2) - - - - - + meanD = (Dpar + 2 * Dperp) / 3. + return np.sqrt(3 / 2.) * np.sqrt((Dpar - meanD) ** 2 + 2 * (Dperp - meanD) ** 2) / np.sqrt( + Dpar ** 2 + 2 * Dperp ** 2) ## Start point for the optimization @@ -242,21 +213,17 @@ def FA(Dpar, Dperp): ## by principle, i'm happier when I initialize my optimization away from my constraint boundaries x0 = np.array([2.0, 1.3, 0.7, 0.65, 0.3, 0.05]) - - - ## storage function_residual = np.zeros(data_SM.shape[:3]) exit_code = np.zeros(data_SM.shape[:3]) -fitted_param = np.zeros(data_SM.shape[:3]+(6,)) ## this would be (5,) for model1 - +fitted_param = np.zeros(data_SM.shape[:3] + (6,)) ## this would be (5,) for model1 # voxelwise minimization with linear and ratio constraint def run_optim_fixed_x0(datavox): - res = opt.minimize(func_vox_ls, x0=x0, args=(b, bd, datavox), jac='3-point', hess=opt.BFGS(), method='trust-constr', constraints=[cons, fa_cons]) - return res - + res = opt.minimize(func_vox_ls, x0=x0, args=(b, bd, datavox), jac='3-point', hess=opt.BFGS(), method='trust-constr', + constraints=[cons, fa_cons]) + return res ## loop over voxels in the mask with multiprocess.Pool.map @@ -272,44 +239,34 @@ def run_optim_fixed_x0(datavox): print('time = {} s\n'.format(t2 - t1)) - ## store values function_residual[mask] = np.array([r.fun for r in res]) exit_code[mask] = np.array([r.status for r in res]) fitted_param[mask] = np.array([r.x for r in res]) - ## save endtag = '1' nib.save(nib.nifti1.Nifti1Image(fitted_param, b0_ima.affine, b0_ima.header), savepath + 'fit{}.nii.gz'.format(endtag)) nib.save(nib.nifti1.Nifti1Image(exit_code, b0_ima.affine, b0_ima.header), savepath + 'exit{}.nii.gz'.format(endtag)) -nib.save(nib.nifti1.Nifti1Image(function_residual, b0_ima.affine, b0_ima.header), savepath + 'res{}.nii.gz'.format(endtag)) - +nib.save(nib.nifti1.Nifti1Image(function_residual, b0_ima.affine, b0_ima.header), + savepath + 'res{}.nii.gz'.format(endtag)) ## compute uFA voxelwise and save uFA = np.zeros(data_SM.shape[:3]) for vox in np.ndindex(data_SM.shape[:3]): - if mask[vox]: - uFA[vox] = microFA(*fitted_param[vox]) + if mask[vox]: + uFA[vox] = microFA(*fitted_param[vox]) nib.save(nib.nifti1.Nifti1Image(uFA, b0_ima.affine, b0_ima.header), savepath + 'uFA{}.nii.gz'.format(endtag)) - - ## compute FA of extra compartement (sanity check / necessity check of nonlin constraint) voxelwise and save FAextra = np.zeros(data_SM.shape[:3]) for vox in np.ndindex(data_SM.shape[:3]): - if mask[vox]: - FAextra[vox] = FA(fitted_param[vox][1],fitted_param[vox][2]) - + if mask[vox]: + FAextra[vox] = FA(fitted_param[vox][1], fitted_param[vox][2]) nib.save(nib.nifti1.Nifti1Image(FAextra, b0_ima.affine, b0_ima.header), savepath + 'FA_extra{}.nii.gz'.format(endtag)) - print('Done') pool.close() - - - - diff --git a/example_scripts/validate_formulas_example.py b/example_scripts/validate_formulas_example.py index f0eab84..9b7c27e 100644 --- a/example_scripts/validate_formulas_example.py +++ b/example_scripts/validate_formulas_example.py @@ -4,7 +4,6 @@ # from sphHist import plotScatter3 from linear_planar.smt_lin_pla import * - ## we want to generate a set of well distribution point on the sphere ## and then generate a linear and a planar B-tensor for each point ## the linear one points in the generated direction @@ -19,87 +18,88 @@ theta = golden_angle * np.arange(N) z = np.linspace(1 - 1.0 / N, 1.0 / N - 1, N) radius = np.sqrt(1 - z * z) - -points = np.zeros((N, 3)) -points[:,0] = radius * np.cos(theta) -points[:,1] = radius * np.sin(theta) -points[:,2] = z +points = np.zeros((N, 3)) +points[:, 0] = radius * np.cos(theta) +points[:, 1] = radius * np.sin(theta) +points[:, 2] = z ## convert to dipy hemisphere hemi = HemiSphere(xyz=points) res, _ = disperse_charges(hemi, 100) dirs = res.vertices - ## get necessary vector to generate the plane planar_dirs = np.zeros((dirs.shape[0], 2, 3)) for idx in range(dirs.shape[0]): - # current directions - v = dirs[idx] - # generate initial non colinear vectors, v1 is v - v1, v2, v3 = findLinIndepRandomRot(v) - # generate 3 orthonormal vectors, u1 is v normalized - u1, u2, u3 = gramSchmidt3(v1, v2, v3) - planar_dirs[idx, 0] = u2 - planar_dirs[idx, 1] = u3 - + # current directions + v = dirs[idx] + # generate initial non colinear vectors, v1 is v + v1, v2, v3 = findLinIndepRandomRot(v) + # generate 3 orthonormal vectors, u1 is v normalized + u1, u2, u3 = gramSchmidt3(v1, v2, v3) + planar_dirs[idx, 0] = u2 + planar_dirs[idx, 1] = u3 ## G1G2G3.txt-like structure ## g3 = 0 because we are not doing spherical B-tensor -g1 = np.concatenate((dirs, planar_dirs[:,0])) -g2 = np.concatenate((dirs, planar_dirs[:,1])) +g1 = np.concatenate((dirs, planar_dirs[:, 0])) +g2 = np.concatenate((dirs, planar_dirs[:, 1])) g3 = np.zeros_like(g1) -grad = np.concatenate((g1,g2,g3), axis=1) +grad = np.concatenate((g1, g2, g3), axis=1) ## 3-gradient form to B-tensor -def table2ten(gtable, bval = 1000.): - tensors_1 = np.zeros((gtable.shape[0], 3, 3)) - tensors_2 = np.zeros((gtable.shape[0], 3, 3)) - tensors_3 = np.zeros((gtable.shape[0], 3, 3)) - - for ix in range(3): - for iy in range(3): - tensors_1[:, ix, iy] = gtable[:, ix]*gtable[:, iy] - tensors_2[:, ix, iy] = gtable[:, ix+3]*gtable[:, iy+3] - tensors_3[:, ix, iy] = gtable[:, ix+6]*gtable[:, iy+6] +def table2ten(gtable, bval=1000.): + tensors_1 = np.zeros((gtable.shape[0], 3, 3)) + tensors_2 = np.zeros((gtable.shape[0], 3, 3)) + tensors_3 = np.zeros((gtable.shape[0], 3, 3)) - ten = bval * (tensors_1 + tensors_2 + tensors_3) - return ten + for ix in range(3): + for iy in range(3): + tensors_1[:, ix, iy] = gtable[:, ix] * gtable[:, iy] + tensors_2[:, ix, iy] = gtable[:, ix + 3] * gtable[:, iy + 3] + tensors_3[:, ix, iy] = gtable[:, ix + 6] * gtable[:, iy + 6] + ten = bval * (tensors_1 + tensors_2 + tensors_3) + return ten ## split shape ## NOTE, the bval is a PER GRADIENT scalling ## so here, we have bval=500 and the resulting true bvalue will be 1000 because g1 and g2 have norm 1 and g3 has norm 0 -ten_lin = table2ten(grad[:N], bval = 500.) -ten_pla = table2ten(grad[N:], bval = 500.) +ten_lin = table2ten(grad[:N], bval=500.) +ten_pla = table2ten(grad[N:], bval=500.) ## convert mm^2/s to um^2/ms -ten_lin = ten_lin/1000. -ten_pla = ten_pla/1000. - - +ten_lin = ten_lin / 1000. +ten_pla = ten_pla / 1000. ## compute bvalues b_lin = np.trace(ten_lin, axis1=1, axis2=2) b_pla = np.trace(ten_pla, axis1=1, axis2=2) - ## necessary tensor-product for signal gen -n = np.array([0,0,1]) # tissue orientation -pb_lin = np.outer(n[0]**2, ten_lin[:,0,0]) + np.outer(n[1]**2, ten_lin[:,1,1]) + np.outer(n[2]**2, ten_lin[:,2,2]) + 2*np.outer(n[0]*n[1], ten_lin[:,0,1]) + 2*np.outer(n[0]*n[2], ten_lin[:,0,2]) + 2*np.outer(n[1]*n[2], ten_lin[:,1,2]) -pb_pla = np.outer(n[0]**2, ten_pla[:,0,0]) + np.outer(n[1]**2, ten_pla[:,1,1]) + np.outer(n[2]**2, ten_pla[:,2,2]) + 2*np.outer(n[0]*n[1], ten_pla[:,0,1]) + 2*np.outer(n[0]*n[2], ten_pla[:,0,2]) + 2*np.outer(n[1]*n[2], ten_pla[:,1,2]) +n = np.array([0, 0, 1]) # tissue orientation +pb_lin = np.outer(n[0] ** 2, ten_lin[:, 0, 0]) + np.outer(n[1] ** 2, ten_lin[:, 1, 1]) + np.outer(n[2] ** 2, + ten_lin[:, 2, + 2]) + 2 * np.outer( + n[0] * n[1], ten_lin[:, 0, 1]) + 2 * np.outer(n[0] * n[2], ten_lin[:, 0, 2]) + 2 * np.outer(n[1] * n[2], + ten_lin[:, 1, 2]) +pb_pla = np.outer(n[0] ** 2, ten_pla[:, 0, 0]) + np.outer(n[1] ** 2, ten_pla[:, 1, 1]) + np.outer(n[2] ** 2, + ten_pla[:, 2, + 2]) + 2 * np.outer( + n[0] * n[1], ten_pla[:, 0, 1]) + 2 * np.outer(n[0] * n[2], ten_pla[:, 0, 2]) + 2 * np.outer(n[1] * n[2], + ten_pla[:, 1, 2]) ## non spherical mean signal generation ## V, Vex, Vw := fractions for intra;extra;freewater def genSig(D1, D2, D3, V, Vex, Vw, pb, b): - S_i = np.exp(-pb*D1[:,None]) - S_e = np.exp(-pb*D2[:,None] -(b[None,:]-pb)*D3[:,None]) - S_w = np.exp(-b[None,:]*3*np.ones((Vw.shape[0],1))) - return V[:,None]*S_i + Vex[:,None]*S_e + Vw[:,None]*S_w + S_i = np.exp(-pb * D1[:, None]) + S_e = np.exp(-pb * D2[:, None] - (b[None, :] - pb) * D3[:, None]) + S_w = np.exp(-b[None, :] * 3 * np.ones((Vw.shape[0], 1))) + return V[:, None] * S_i + Vex[:, None] * S_e + Vw[:, None] * S_w ## test values @@ -110,27 +110,20 @@ def genSig(D1, D2, D3, V, Vex, Vw, pb, b): Vex = np.array([1]) Vw = np.array([0]) - S_lin = genSig(D1, D2, D3, V, Vex, Vw, pb_lin, b_lin) S_pla = genSig(D1, D2, D3, V, Vex, Vw, pb_pla, b_pla) - ## this should match signal_lin_numerical = S_lin.mean() signal_lin_analytical = sm_signal_model_2(b_lin[0], 1., D1, D2, D3, V, Vex, Vw)[0] signal_pla_numerical = S_pla.mean() signal_pla_analytical = sm_signal_model_2(b_pla[0], -0.5, D1, D2, D3, V, Vex, Vw)[0] - print('\n') print('signal linear numerical {}'.format(signal_lin_numerical)) print('signal linear analytical {}'.format(signal_lin_analytical)) -print('abs diff {:.2e}'.format(np.abs(signal_lin_numerical-signal_lin_analytical))) +print('abs diff {:.2e}'.format(np.abs(signal_lin_numerical - signal_lin_analytical))) print('\n') print('signal linear numerical {}'.format(signal_pla_numerical)) print('signal linear analytical {}'.format(signal_pla_analytical)) -print('abs diff {:.2e}'.format(np.abs(signal_pla_numerical-signal_pla_analytical))) - - - - +print('abs diff {:.2e}'.format(np.abs(signal_pla_numerical - signal_pla_analytical))) diff --git a/linear_planar/calc_jacobian.py b/linear_planar/calc_jacobian.py index fecc957..8cb4504 100644 --- a/linear_planar/calc_jacobian.py +++ b/linear_planar/calc_jacobian.py @@ -10,37 +10,37 @@ # Hard-coded max jacobian determinant for siemens, taken from grad unwarp file of Human Connectom Project siemens_max_det = 10. -DESCRIPTION = 'Calculation of the Jacobi Determinant from a given warp field. Cornelius Eichner 2018' +DESCRIPTION = 'Calculation of the Jacobi Determinant from a given warp field. Cornelius Eichner 2018' + def buildArgsParser(): p = argparse.ArgumentParser(description=DESCRIPTION) p.add_argument('--in', dest='data', action='store', type=str, - help='Path of the input warp field (nifti format)') + help='Path of the input warp field (nifti format)') p.add_argument('--out', dest='out', action='store', type=str, - help='Path of the output') + help='Path of the output') return p def evaluate_jacobian(F, pixdim): - # Differentials in each dimension d0, d1, d2 = pixdim[0], pixdim[1], pixdim[2] if d0 == 0 or d1 == 0 or d2 == 0: raise ValueError('weirdness found in Jacobian calculation') - dFxdx, dFxdy, dFxdz = np.gradient(F[...,0], d0, d1, d2) - dFydx, dFydy, dFydz = np.gradient(F[...,1], d0, d1, d2) - dFzdx, dFzdy, dFzdz = np.gradient(F[...,2], d0, d1, d2) + dFxdx, dFxdy, dFxdz = np.gradient(F[..., 0], d0, d1, d2) + dFydx, dFydy, dFydz = np.gradient(F[..., 1], d0, d1, d2) + dFzdx, dFzdy, dFzdz = np.gradient(F[..., 2], d0, d1, d2) jacdet = (dFxdx) * (dFydy) * (dFzdz) \ - - (dFxdx) * dFydz * dFzdy \ - - dFxdy * dFydx * (dFzdz) \ - + dFxdy * dFydz * dFzdx \ - + dFxdz * dFydx * dFzdy \ - - dFxdz * (dFydy) * dFzdx + - (dFxdx) * dFydz * dFzdy \ + - dFxdy * dFydx * (dFzdz) \ + + dFxdy * dFydz * dFzdx \ + + dFxdz * dFydx * dFzdy \ + - dFxdz * (dFydy) * dFzdx jacdet = np.abs(jacdet) jacdet[np.where(jacdet > siemens_max_det)] = siemens_max_det @@ -48,7 +48,6 @@ def evaluate_jacobian(F, pixdim): return jacdet - def main(): # Load parser to read data from command line input parser = buildArgsParser() diff --git a/linear_planar/eddy_nlgc_mp.py b/linear_planar/eddy_nlgc_mp.py index d919f12..a84dbbe 100644 --- a/linear_planar/eddy_nlgc_mp.py +++ b/linear_planar/eddy_nlgc_mp.py @@ -6,7 +6,7 @@ import nibabel as nib import numpy as np -DESCRIPTION = 'Joint correction of eddy currents and gradient non linearity correction using FSL tools. Cornelius Eichner 2018' +DESCRIPTION = 'Joint correction of eddy currents and gradient non linearity correction using FSL tools. Cornelius Eichner 2018' PATH_GRAD_UNWARP = '/home/raid2/paquette/tools/gradunwarp/gradunwarp/core/gradient_unwarp.py' # PATH_EDDY = '/home/raid2/paquette/tools/eddy_openmp' @@ -17,37 +17,37 @@ def buildArgsParser(): p = argparse.ArgumentParser(description=DESCRIPTION) p.add_argument('--in', dest='data', action='store', type=str, - help='Path of the input volume (nifti format)') + help='Path of the input volume (nifti format)') p.add_argument('--bvec', dest='bvec', action='store', type=str, - help='Path of the bvec file (default \'bvec\')') + help='Path of the bvec file (default \'bvec\')') p.add_argument('--bval', dest='bval', action='store', type=str, - help='Path of the bval file (default \'bval\')') + help='Path of the bval file (default \'bval\')') p.add_argument('--mask', dest='mask', action='store', type=str, - help='Path of the brain mask') + help='Path of the brain mask') p.add_argument('--acqp', dest='acqp', action='store', type=str, - help='Path of eddy acquisition parameter file (default \'acq_param.txt\')') - + help='Path of eddy acquisition parameter file (default \'acq_param.txt\')') + p.add_argument('--index', dest='index', action='store', type=str, - help='Path of the eddy index file (default \'index.txt\')') + help='Path of the eddy index file (default \'index.txt\')') p.add_argument('--topup', dest='topup', action='store', type=str, - help='Base name of topup file structure (default \'topup/topup\')') + help='Base name of topup file structure (default \'topup/topup\')') p.add_argument('--out', dest='out', action='store', type=str, - help='Path of the output volume') + help='Path of the output volume') p.add_argument('--pmat', dest='pmat', action='store', type=str, - help='Path of the affine transform to apply at the end of registrationb (e.g., T1 anatomy)') + help='Path of the affine transform to apply at the end of registrationb (e.g., T1 anatomy)') - p.add_argument('--mp', dest='openmp', action='store', type=int, default='1', - help='Optional: OpenMP parallelization factor (default 8)') + p.add_argument('--mp', dest='openmp', action='store', type=int, default='1', + help='Optional: OpenMP parallelization factor (default 8)') p.add_argument('--jump', dest='jump', action='store', type=str, - help='skip part of the processing') + help='skip part of the processing') return p @@ -107,9 +107,8 @@ def main(): DATA = os.path.realpath(args.data) MASK = os.path.realpath(args.mask) - - OUT = args.out + OUT = args.out ###################### # Calculate Gradient Non Linearities @@ -121,9 +120,9 @@ def main(): os.system(cmd) # Perform gradient non linearity correction - cmd = PATH_GRAD_UNWARP + ' '+ PATH + 'single.nii.gz ' + PATH + 'nlgc.nii.gz siemens -g /home/raid2/paquette/tools/gradunwarp/coeffConnectom.grad -n' + cmd = PATH_GRAD_UNWARP + ' ' + PATH + 'single.nii.gz ' + PATH + 'nlgc.nii.gz siemens -g /home/raid2/paquette/tools/gradunwarp/coeffConnectom.grad -n' os.system(cmd) - + # Rename resulting files cmd = 'mv fullWarp_abs.nii.gz ' + PATH + 'nlgc_warp.nii.gz' os.system(cmd) @@ -137,7 +136,6 @@ def main(): nlgc_warp.nii.gz warp field for gradient correction """ - ###################### # Calculate eddy displacement fields if jump < 2: @@ -148,9 +146,9 @@ def main(): os.system(cmd) # Run eddy from home installation - cmd =\ - PATH_EDDY + ' \ - --imain='+ DATA + ' \ + cmd = \ + PATH_EDDY + ' \ + --imain=' + DATA + ' \ --mask=' + MASK + ' \ --index=' + INDEX + ' \ --acqp=' + ACQP + ' \ @@ -166,7 +164,6 @@ def main(): os.system(cmd) os.system("echo Eddy Complete") - # Move all displacement field files in dfield subdirectory os.system("echo Moving displacement field files into common folder") os.system('mkdir -p dfield') @@ -185,7 +182,6 @@ def main(): dfield/ Folder containing all final deformation fields from eddy """ - ###################### # Combine both warp fields if jump < 3: @@ -221,7 +217,6 @@ def main(): comb_warp/ Folder containing all combined warp files """ - ###################### # Apply warp for using combined field os.system("echo Application of combined warp fields") @@ -248,7 +243,6 @@ def main(): else: print('jumping over applywarp') - ###################### # Correct for warp signal bias using Jacobian determinante os.system("echo Correction of signal intensities with Jacobian Determinant") @@ -257,13 +251,14 @@ def main(): os.system('mkdir -p corr_data') for i in xrange(len(os.listdir('comb_warp/'))): - cmd = 'python ' + PATH_CALC_JACOBIAN + ' --in comb_warp/' + COMB_WARP_FILES[i] + ' --out jac_det_data/jac.' + COMB_WARP_FILES[i] + cmd = 'python ' + PATH_CALC_JACOBIAN + ' --in comb_warp/' + COMB_WARP_FILES[i] + ' --out jac_det_data/jac.' + \ + COMB_WARP_FILES[i] os.system(cmd) - cmd = 'fslmaths warped_data/' + SPLIT_DATA_FILES[i] + ' -mul jac_det_data/jac.' + COMB_WARP_FILES[i] + ' corr_data/' + SPLIT_DATA_FILES[i] + cmd = 'fslmaths warped_data/' + SPLIT_DATA_FILES[i] + ' -mul jac_det_data/jac.' + COMB_WARP_FILES[ + i] + ' corr_data/' + SPLIT_DATA_FILES[i] # cmd = 'fslmaths warped_data/' + SPLIT_DATA_FILES[i] + ' -div jac_det_data/jac.' + COMB_WARP_FILES[i] + ' corr_data/' + SPLIT_DATA_FILES[i] os.system(cmd) - ###################### # Merge warped volumes os.system("echo Merging corrected volumes") diff --git a/linear_planar/gs.py b/linear_planar/gs.py index b562b89..fd3ad53 100644 --- a/linear_planar/gs.py +++ b/linear_planar/gs.py @@ -1,59 +1,62 @@ import numpy as np -def proj(u,v): - return (np.dot(u,v)/np.dot(u,u)) * u + +def proj(u, v): + return (np.dot(u, v) / np.dot(u, u)) * u + def angleVec3(v1, v2): - ang = np.math.acos(np.clip(np.dot(v1,v2)/(np.linalg.norm(v1)*np.linalg.norm(v2)),-1,1)) - return ang + ang = np.math.acos(np.clip(np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2)), -1, 1)) + return ang + def checkOrtho3(u1, u2, u3, tol=1e-3): - a12 = angleVec3(u1, u2) - a13 = angleVec3(u1, u3) - a23 = angleVec3(u2, u3) - if (np.abs(a12-np.pi/2.) < tol) and (np.abs(a13-np.pi/2.) < tol) and (np.abs(a23-np.pi/2.) < tol): - return True - else: - return False + a12 = angleVec3(u1, u2) + a13 = angleVec3(u1, u3) + a23 = angleVec3(u2, u3) + if (np.abs(a12 - np.pi / 2.) < tol) and (np.abs(a13 - np.pi / 2.) < tol) and (np.abs(a23 - np.pi / 2.) < tol): + return True + else: + return False + def findLinIndep(v): - # Assumes the input has norm > 0 - v1 = v / np.linalg.norm(v) - idx = np.argsort(v1) - v2 = np.zeros(3) - v2[idx[0]] = 1 - v3 = np.zeros(3) - v3[idx[1]] = 1 - return v1, v2, v3 + # Assumes the input has norm > 0 + v1 = v / np.linalg.norm(v) + idx = np.argsort(v1) + v2 = np.zeros(3) + v2[idx[0]] = 1 + v3 = np.zeros(3) + v3[idx[1]] = 1 + return v1, v2, v3 def findLinIndepRandomRot(v): - # Assumes the input has norm > 0 - v1 = v / np.linalg.norm(v) - idx = np.argsort(v1) - ang = 2*np.pi*np.random.rand() - v2 = np.zeros(3) - v2[idx[0]] = np.cos(ang) - v2[idx[1]] = np.sin(ang) - v3 = np.zeros(3) - v3[idx[0]] = -np.sin(ang) - v3[idx[1]] = np.cos(ang) - return v1, v2, v3 + # Assumes the input has norm > 0 + v1 = v / np.linalg.norm(v) + idx = np.argsort(v1) + ang = 2 * np.pi * np.random.rand() + v2 = np.zeros(3) + v2[idx[0]] = np.cos(ang) + v2[idx[1]] = np.sin(ang) + v3 = np.zeros(3) + v3[idx[0]] = -np.sin(ang) + v3[idx[1]] = np.cos(ang) + return v1, v2, v3 def gramSchmidt3(v1, v2, v3): - # produce vectors u1 u2 u3 which are orthogonal and e1 e2 e3 which are orthonormal - # all vector are normalized - # need v1 v2 v3 to be linearly independant (use findLinIndep()) - # u1 has the same orienation as v1 - # Gram-Schmidt 3 vector - u1 = v1.copy() - u2 = v2 - proj(u1,v2) - u3 = v3 - proj(u1,v3) - proj(u2,v3) - # e1,e2,e3 are orthonormal - # e1 is parallel to dir therefore e2,e3 span a perpendicular plane to dir - e1 = u1 / np.linalg.norm(u1) - e2 = u2 / np.linalg.norm(u2) - e3 = u3 / np.linalg.norm(u3) - return e1, e2, e3 - + # produce vectors u1 u2 u3 which are orthogonal and e1 e2 e3 which are orthonormal + # all vector are normalized + # need v1 v2 v3 to be linearly independant (use findLinIndep()) + # u1 has the same orienation as v1 + # Gram-Schmidt 3 vector + u1 = v1.copy() + u2 = v2 - proj(u1, v2) + u3 = v3 - proj(u1, v3) - proj(u2, v3) + # e1,e2,e3 are orthonormal + # e1 is parallel to dir therefore e2,e3 span a perpendicular plane to dir + e1 = u1 / np.linalg.norm(u1) + e2 = u2 / np.linalg.norm(u2) + e3 = u3 / np.linalg.norm(u3) + return e1, e2, e3 diff --git a/linear_planar/mean_std_4th.py b/linear_planar/mean_std_4th.py index bc542aa..8c11db2 100644 --- a/linear_planar/mean_std_4th.py +++ b/linear_planar/mean_std_4th.py @@ -3,19 +3,19 @@ import sys if __name__ == "__main__": - print('\nparam order::\ninput_vol.nii(.gz) output_vol_mean.nii(.gz) output_vol_std.nii(.gz)\n') + print('\nparam order::\ninput_vol.nii(.gz) output_vol_mean.nii(.gz) output_vol_std.nii(.gz)\n') - # loading data - img = nib.load(sys.argv[1]) - data = img.get_data() + # loading data + img = nib.load(sys.argv[1]) + data = img.get_data() - newdata1 = data.mean(axis=3) - newdata2 = data.std(axis=3) + newdata1 = data.mean(axis=3) + newdata2 = data.std(axis=3) - newimg1 = nib.nifti1.Nifti1Image(newdata1, img.affine, img.header) - newimg2 = nib.nifti1.Nifti1Image(newdata2, img.affine, img.header) + newimg1 = nib.nifti1.Nifti1Image(newdata1, img.affine, img.header) + newimg2 = nib.nifti1.Nifti1Image(newdata2, img.affine, img.header) - nib.save(newimg1, sys.argv[2]) - nib.save(newimg2, sys.argv[3]) + nib.save(newimg1, sys.argv[2]) + nib.save(newimg2, sys.argv[3]) diff --git a/linear_planar/smt_lin_pla.py b/linear_planar/smt_lin_pla.py index 7dd2ba9..6c26dc0 100644 --- a/linear_planar/smt_lin_pla.py +++ b/linear_planar/smt_lin_pla.py @@ -1,4 +1,4 @@ -import numpy as np +import numpy as np from scipy.special import erf import scipy as sp @@ -39,70 +39,70 @@ # b_delta = bd = {-0.5 if planar; 0 if spherical; 1 if linear} - def main(): - # a few b_delta - bd_pla = -0.5 - bd_sph = 0.0 - bd_lin = 1.0 - - return None + # a few b_delta + bd_pla = -0.5 + bd_sph = 0.0 + bd_lin = 1.0 + return None # param conversion utils def Di_from_param(Dpar, Dperp): - return (Dpar+2*Dperp)/3. + return (Dpar + 2 * Dperp) / 3. + def Dd_from_param(Dpar, Dperp): - return (Dpar-Dperp)/float(Dpar+2*Dperp) + return (Dpar - Dperp) / float(Dpar + 2 * Dperp) -def Dpar_from_param(Di, Dd): - return Di*(1+2*Dd) -def Dperp_from_param(Di, Dd): - return Di*(1-Dd) +def Dpar_from_param(Di, Dd): + return Di * (1 + 2 * Dd) +def Dperp_from_param(Di, Dd): + return Di * (1 - Dd) def gfunc(alpha): - # return sp.sqrt(np.pi/(4*alpha)) * erf(sp.sqrt(alpha)) - tmp = sp.sqrt(np.pi/np.float64(4.*alpha)) * erf(sp.sqrt(alpha)) - return np.where(np.abs(alpha)<1e-16, 1., np.abs(tmp)) # the abs is missing from the paper + # return sp.sqrt(np.pi/(4*alpha)) * erf(sp.sqrt(alpha)) + tmp = sp.sqrt(np.pi / np.float64(4. * alpha)) * erf(sp.sqrt(alpha)) + return np.where(np.abs(alpha) < 1e-16, 1., np.abs(tmp)) # the abs is missing from the paper + # generic signal formula for any B-tensor shape and any D-tensor shape def sm_signal_generic(b, bd, Di, Dd): - return np.exp(-b*Di*(1-bd*Dd)) * gfunc(3*b*Di*bd*Dd) + return np.exp(-b * Di * (1 - bd * Dd)) * gfunc(3 * b * Di * bd * Dd) + # generic signal formula for any B-tensor shape and any D-tensor shape reparametrized def sm_signal_generic_reparam(b, bd, Dpar, Dperp): - Di = Di_from_param(Dpar, Dperp) - Dd = Dd_from_param(Dpar, Dperp) - return sm_signal_generic(b, bd, Di, Dd) + Di = Di_from_param(Dpar, Dperp) + Dd = Dd_from_param(Dpar, Dperp) + return sm_signal_generic(b, bd, Di, Dd) # model1 stick(intra) + zepplin(extra) + ball(csf) def sm_signal_model_1(b, bd, Dpar_in, Dpar_ex, Dperp_ex, f_in, f_ex): - # stick-like intra axonal compartement with Dperp_ex = 0 - sm_signal_in = sm_signal_generic_reparam(b, bd, Dpar_in, 0) - # zepplin-like extra axonal compartement - sm_signal_ex = sm_signal_generic_reparam(b, bd, Dpar_ex, Dperp_ex) - # ball-like csf axonal compartement with Dpar_ex = Dperp_ex = 3 um^2/ms (free water diffusivity at in-vivo human brain temperature) - sm_signal_csf = sm_signal_generic_reparam(b, bd, 3., 3.) - # model1 enforces sum(fraction) = 1 indirectly - f_csf = 1 - f_in - f_ex - return f_in*sm_signal_in + f_ex*sm_signal_ex + f_csf*sm_signal_csf + # stick-like intra axonal compartement with Dperp_ex = 0 + sm_signal_in = sm_signal_generic_reparam(b, bd, Dpar_in, 0) + # zepplin-like extra axonal compartement + sm_signal_ex = sm_signal_generic_reparam(b, bd, Dpar_ex, Dperp_ex) + # ball-like csf axonal compartement with Dpar_ex = Dperp_ex = 3 um^2/ms (free water diffusivity at in-vivo human brain temperature) + sm_signal_csf = sm_signal_generic_reparam(b, bd, 3., 3.) + # model1 enforces sum(fraction) = 1 indirectly + f_csf = 1 - f_in - f_ex + return f_in * sm_signal_in + f_ex * sm_signal_ex + f_csf * sm_signal_csf + # model2 stick(intra) + zepplin(extra) + ball(csf) def sm_signal_model_2(b, bd, Dpar_in, Dpar_ex, Dperp_ex, f_in, f_ex, f_csf): - # stick-like intra axonal compartement with Dperp_ex = 0 - sm_signal_in = sm_signal_generic_reparam(b, bd, Dpar_in, 0) - # zepplin-like extra axonal compartement - sm_signal_ex = sm_signal_generic_reparam(b, bd, Dpar_ex, Dperp_ex) - # ball-like csf axonal compartement with Dpar_ex = Dperp_ex = 3 um^2/ms (free water diffusivity at in-vivo human brain temperature) - sm_signal_csf = sm_signal_generic_reparam(b, bd, 3., 3.) - # model2 doesnt enforces sum(fraction) = 1 directly, need to deal with it "later" - return f_in*sm_signal_in + f_ex*sm_signal_ex + f_csf*sm_signal_csf - - + # stick-like intra axonal compartement with Dperp_ex = 0 + sm_signal_in = sm_signal_generic_reparam(b, bd, Dpar_in, 0) + # zepplin-like extra axonal compartement + sm_signal_ex = sm_signal_generic_reparam(b, bd, Dpar_ex, Dperp_ex) + # ball-like csf axonal compartement with Dpar_ex = Dperp_ex = 3 um^2/ms (free water diffusivity at in-vivo human brain temperature) + sm_signal_csf = sm_signal_generic_reparam(b, bd, 3., 3.) + # model2 doesnt enforces sum(fraction) = 1 directly, need to deal with it "later" + return f_in * sm_signal_in + f_ex * sm_signal_ex + f_csf * sm_signal_csf diff --git a/linear_planar/split_shells_and_types.py b/linear_planar/split_shells_and_types.py index 7c9a118..5246cef 100644 --- a/linear_planar/split_shells_and_types.py +++ b/linear_planar/split_shells_and_types.py @@ -3,65 +3,66 @@ import sys -def main(dwifile, bvals, gtablefile, bmax, eddy_configfile, eddy_indicesfile, dwiout, bvecout, bvalout, eddy_indices_out, position_idx_out): +def main(dwifile, bvals, gtablefile, bmax, eddy_configfile, eddy_indicesfile, dwiout, bvecout, bvalout, + eddy_indices_out, position_idx_out): + # convert cmdline strings + bvals = [int(i) for i in bvals.split(',')] + bmax = int(bmax) - # convert cmdline strings - bvals = [int(i) for i in bvals.split(',')] - bmax = int(bmax) + # load data + dwiimg = nib.load(dwifile) + dwiall = dwiimg.get_data() - # load data - dwiimg = nib.load(dwifile) - dwiall = dwiimg.get_data() + gtable = np.genfromtxt(gtablefile) + eddy_indices = np.genfromtxt(eddy_indicesfile) - gtable = np.genfromtxt(gtablefile) - eddy_indices = np.genfromtxt(eddy_indicesfile) + # guess bval and type from gradient table + gnorm = np.linalg.norm(gtable[:, :3], axis=1) + # assumes g1 and g2 have same bvalue + databval = bmax * gnorm ** 2 - # guess bval and type from gradient table - gnorm = np.linalg.norm(gtable[:, :3], axis=1) - # assumes g1 and g2 have same bvalue - databval = bmax*gnorm**2 + # linear-planar detection (assumes only perfect linear and planar exist (eucl. dist. between g1 and g2 is 0 OR + # sqrt(2)~1.4)) + isPlanar = np.linalg.norm((gtable[:, :3] - gtable[:, 3:6]) / gnorm[:, None], axis=1) > 1 + # detect which shell in bvals fits best dataval + shell_idx = np.argmin(np.abs(databval - np.array(bvals)[:, None]), axis=0) - # linear-planar detection (assumes only perfect linear and planar exist (eucl. dist. between g1 and g2 is 0 OR sqrt(2)~1.4)) - isPlanar = np.linalg.norm((gtable[:, :3]-gtable[:, 3:6])/gnorm[:,None], axis=1) > 1 - # detect which shell in bvals fits best dataval - shell_idx = np.argmin(np.abs(databval - np.array(bvals)[:,None]), axis=0) + # split it + for idx, bval in enumerate(bvals): + # checks for linear, if any + b_mask = np.logical_and(shell_idx == idx, isPlanar == 0) + if np.sum(b_mask) > 0: + print('{} linear b = {} volumes'.format(np.sum(b_mask), bval)) + position_idx = np.where(b_mask)[0] + np.savetxt(position_idx_out + '__b{}_lin.txt'.format(bval), position_idx, fmt='%i') + np.savetxt(eddy_indices_out + '__b{}_lin.txt'.format(bval), eddy_indices[position_idx], fmt='%i') + np.savetxt(bvalout + '__b{}_lin.txt'.format(bval), databval[position_idx]) + dirs = gtable[position_idx, :3] + normdirs = dirs / np.linalg.norm(dirs, axis=1)[:, None] + np.savetxt(bvecout + '__b{}_lin.txt'.format(bval), normdirs) + newimg = nib.nifti1.Nifti1Image(dwiall[..., position_idx], dwiimg.affine, dwiimg.header) + nib.save(newimg, dwiout + '__b{}_lin.nii'.format(bval)) + newimg2 = nib.nifti1.Nifti1Image(dwiall[..., position_idx].mean(3), dwiimg.affine, dwiimg.header) + nib.save(newimg2, dwiout + '__b{}_lin_SM.nii'.format(bval)) - # split it - for idx, bval in enumerate(bvals): - # checks for linear, if any - b_mask = np.logical_and(shell_idx==idx, isPlanar==0) - if np.sum(b_mask) > 0: - print('{} linear b = {} volumes'.format(np.sum(b_mask), bval)) - position_idx = np.where(b_mask)[0] - np.savetxt(position_idx_out + '__b{}_lin.txt'.format(bval), position_idx, fmt='%i') - np.savetxt(eddy_indices_out + '__b{}_lin.txt'.format(bval), eddy_indices[position_idx], fmt='%i') - np.savetxt(bvalout + '__b{}_lin.txt'.format(bval), databval[position_idx]) - dirs = gtable[position_idx, :3] - normdirs = dirs / np.linalg.norm(dirs, axis=1)[:,None] - np.savetxt(bvecout + '__b{}_lin.txt'.format(bval), normdirs) - newimg = nib.nifti1.Nifti1Image(dwiall[...,position_idx], dwiimg.affine, dwiimg.header) - nib.save(newimg, dwiout + '__b{}_lin.nii'.format(bval)) - newimg2 = nib.nifti1.Nifti1Image(dwiall[...,position_idx].mean(3), dwiimg.affine, dwiimg.header) - nib.save(newimg2, dwiout + '__b{}_lin_SM.nii'.format(bval)) + # checks for planar, if any + b_mask = np.logical_and(shell_idx == idx, isPlanar == 1) + if np.sum(b_mask) > 0: + print('{} planar b = {} volumes'.format(np.sum(b_mask), bval)) + position_idx = np.where(b_mask)[0] + np.savetxt(position_idx_out + '__b{}_pla.txt'.format(bval), position_idx, fmt='%i') + np.savetxt(eddy_indices_out + '__b{}_pla.txt'.format(bval), eddy_indices[position_idx], fmt='%i') + np.savetxt(bvalout + '__b{}_pla.txt'.format(bval), databval[position_idx]) + dirs = gtable[position_idx, :3] + gtable[position_idx, 3:6] + normdirs = dirs / np.linalg.norm(dirs, axis=1)[:, None] + np.savetxt(bvecout + '__b{}_pla.txt'.format(bval), normdirs) + newimg = nib.nifti1.Nifti1Image(dwiall[..., position_idx], dwiimg.affine, dwiimg.header) + nib.save(newimg, dwiout + '__b{}_pla.nii'.format(bval)) + newimg2 = nib.nifti1.Nifti1Image(dwiall[..., position_idx].mean(3), dwiimg.affine, dwiimg.header) + nib.save(newimg2, dwiout + '__b{}_pla_SM.nii'.format(bval)) - # checks for planar, if any - b_mask = np.logical_and(shell_idx==idx, isPlanar==1) - if np.sum(b_mask) > 0: - print('{} planar b = {} volumes'.format(np.sum(b_mask), bval)) - position_idx = np.where(b_mask)[0] - np.savetxt(position_idx_out + '__b{}_pla.txt'.format(bval), position_idx, fmt='%i') - np.savetxt(eddy_indices_out + '__b{}_pla.txt'.format(bval), eddy_indices[position_idx], fmt='%i') - np.savetxt(bvalout + '__b{}_pla.txt'.format(bval), databval[position_idx]) - dirs = gtable[position_idx, :3] + gtable[position_idx, 3:6] - normdirs = dirs / np.linalg.norm(dirs, axis=1)[:,None] - np.savetxt(bvecout + '__b{}_pla.txt'.format(bval), normdirs) - newimg = nib.nifti1.Nifti1Image(dwiall[...,position_idx], dwiimg.affine, dwiimg.header) - nib.save(newimg, dwiout + '__b{}_pla.nii'.format(bval)) - newimg2 = nib.nifti1.Nifti1Image(dwiall[...,position_idx].mean(3), dwiimg.affine, dwiimg.header) - nib.save(newimg2, dwiout + '__b{}_pla_SM.nii'.format(bval)) if __name__ == "__main__": - # for idx,i in enumerate(sys.argv): - # print(idx, i) - main(*sys.argv[1:]) - + # for idx,i in enumerate(sys.argv): + # print(idx, i) + main(*sys.argv[1:]) From 8adc5fbe5d854b021fa78dff65fc504481dc41f3 Mon Sep 17 00:00:00 2001 From: Daniel Rose Date: Mon, 12 Oct 2020 11:44:08 +0200 Subject: [PATCH 3/5] python 3 compatibility and cosmetics --- README.md | 2 ++ linear_planar/debias.py | 2 +- linear_planar/eddy_nlgc_mp.py | 10 ++++++---- linear_planar/smt_lin_pla.py | 12 ++++++------ 4 files changed, 15 insertions(+), 11 deletions(-) diff --git a/README.md b/README.md index dc01a05..5456893 100644 --- a/README.md +++ b/README.md @@ -30,5 +30,7 @@ This repository contains the preprocessing and fitting scripts for double PGSE d +### Steps to make this work: +1. download and install `autodmri` [github](https://github.com/samuelstjean/autodmri) diff --git a/linear_planar/debias.py b/linear_planar/debias.py index ff75067..902b334 100644 --- a/linear_planar/debias.py +++ b/linear_planar/debias.py @@ -9,7 +9,7 @@ def main(fdata, fNs, fsigmas, fcorr_data, fnan_map, fnan_mean): Ns = nib.load(fNs).get_data() sigmas = nib.load(fsigmas).get_data() - corr_data = np.sqrt(data**2 - 2*Ns[...,None]*(sigmas[...,None]**2)) + corr_data = np.sqrt(data ** 2 - 2 * Ns[..., None] * (sigmas[..., None] ** 2)) nan_map = np.zeros_like(corr_data, dtype=np.bool) nan_idx = np.where(np.isnan(corr_data)) diff --git a/linear_planar/eddy_nlgc_mp.py b/linear_planar/eddy_nlgc_mp.py index a84dbbe..7efb351 100644 --- a/linear_planar/eddy_nlgc_mp.py +++ b/linear_planar/eddy_nlgc_mp.py @@ -6,7 +6,8 @@ import nibabel as nib import numpy as np -DESCRIPTION = 'Joint correction of eddy currents and gradient non linearity correction using FSL tools. Cornelius Eichner 2018' +DESCRIPTION = 'Joint correction of eddy currents and gradient non linearity correction using FSL tools. Cornelius ' \ + 'Eichner 2018 ' PATH_GRAD_UNWARP = '/home/raid2/paquette/tools/gradunwarp/gradunwarp/core/gradient_unwarp.py' # PATH_EDDY = '/home/raid2/paquette/tools/eddy_openmp' @@ -231,7 +232,7 @@ def main(): SPLIT_DATA_FILES = sorted(os.listdir('split_data/')) if jump < 4: - for i in xrange(len(os.listdir('split_data/'))): + for i in range(len(os.listdir('split_data/'))): cmd = 'applywarp \ -i split_data/' + SPLIT_DATA_FILES[i] + ' \ -r split_data/' + SPLIT_DATA_FILES[0] + ' \ @@ -250,13 +251,14 @@ def main(): os.system('mkdir -p jac_det_data') os.system('mkdir -p corr_data') - for i in xrange(len(os.listdir('comb_warp/'))): + for i in range(len(os.listdir('comb_warp/'))): cmd = 'python ' + PATH_CALC_JACOBIAN + ' --in comb_warp/' + COMB_WARP_FILES[i] + ' --out jac_det_data/jac.' + \ COMB_WARP_FILES[i] os.system(cmd) cmd = 'fslmaths warped_data/' + SPLIT_DATA_FILES[i] + ' -mul jac_det_data/jac.' + COMB_WARP_FILES[ i] + ' corr_data/' + SPLIT_DATA_FILES[i] - # cmd = 'fslmaths warped_data/' + SPLIT_DATA_FILES[i] + ' -div jac_det_data/jac.' + COMB_WARP_FILES[i] + ' corr_data/' + SPLIT_DATA_FILES[i] + # cmd = 'fslmaths warped_data/' + SPLIT_DATA_FILES[i] + ' -div jac_det_data/jac.' + COMB_WARP_FILES[i] + ' + # corr_data/' + SPLIT_DATA_FILES[i] os.system(cmd) ###################### diff --git a/linear_planar/smt_lin_pla.py b/linear_planar/smt_lin_pla.py index 6c26dc0..ae080db 100644 --- a/linear_planar/smt_lin_pla.py +++ b/linear_planar/smt_lin_pla.py @@ -18,17 +18,17 @@ # A_T2;k is the T2 attenuation of compartement k # A_D;k is the diffusion attenuation of compartement k -# in practice, because we normalize the measured signal by B0 -# you either assume equal T1, T2, PD for each compartement and fit "volume fractions" -# or you "do" something like f'_k = f_k * S_PD;k * A_T1;k * A_T2;k -# where f'_k is a signal fraction instead of volume fraction, (i.e. volume fractions weighted by "some exponential" that depend on the underlying tissues parameters) +# in practice, because we normalize the measured signal by B0 you either assume equal T1, T2, PD for each +# compartement and fit "volume fractions" or you "do" something like f'_k = f_k * S_PD;k * A_T1;k * A_T2;k where f'_k +# is a signal fraction instead of volume fraction, (i.e. volume fractions weighted by "some exponential" that depend +# on the underlying tissues parameters) # Everything here assumes axi-symmetric B-tensor and D-tensor # so for exemple B = diag(b_para b_perp b_perp) (if this specific B-tensor is "pointing" to [1,0,0]) # so for exemple D = diag(d_para d_perp d_perp) (if this specific D-tensor is "pointing" to [1,0,0]) -# because orientation is meaningless for Spherical mean, -# both B-tensor and D-tensor are reparametrize into a "size" parameter and a "shape" parameters (can also be seen as an Isotropic and an Anisotropic parameter) +# because orientation is meaningless for Spherical mean, both B-tensor and D-tensor are reparametrize into a "size" +# parameter and a "shape" parameters (can also be seen as an Isotropic and an Anisotropic parameter) # D_Iso = Di = (D_para + 2*D_perp) / 3 From 0ed13b25c12f30d4305381c23478c69784cdda33 Mon Sep 17 00:00:00 2001 From: Daniel Rose Date: Tue, 10 Nov 2020 13:56:21 +0100 Subject: [PATCH 4/5] refactored command line script into command line tool with function --- linear_planar/mean_std_4th.py | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/linear_planar/mean_std_4th.py b/linear_planar/mean_std_4th.py index 8c11db2..1719fa9 100644 --- a/linear_planar/mean_std_4th.py +++ b/linear_planar/mean_std_4th.py @@ -2,11 +2,10 @@ import numpy as np import sys -if __name__ == "__main__": - print('\nparam order::\ninput_vol.nii(.gz) output_vol_mean.nii(.gz) output_vol_std.nii(.gz)\n') +def mean_std_from_timeseries(in_file: str, out_mean_file: str, out_std_file: str): # loading data - img = nib.load(sys.argv[1]) + img = nib.load(in_file) data = img.get_data() newdata1 = data.mean(axis=3) @@ -15,7 +14,13 @@ newimg1 = nib.nifti1.Nifti1Image(newdata1, img.affine, img.header) newimg2 = nib.nifti1.Nifti1Image(newdata2, img.affine, img.header) - nib.save(newimg1, sys.argv[2]) - nib.save(newimg2, sys.argv[3]) + nib.save(newimg1, out_mean_file) + nib.save(newimg2, out_std_file) + + return out_mean_file, out_std_file +if __name__ == "__main__": + print('\nparam order::\ninput_vol.nii(.gz) output_vol_mean.nii(.gz) output_vol_std.nii(.gz)\n') + + mean_std_from_timeseries(*sys.argv) From bffaf0ad09ab79f0aca1be7541c6cba1fb5ac3cf Mon Sep 17 00:00:00 2001 From: Daniel Rose Date: Tue, 24 Nov 2020 11:38:43 +0100 Subject: [PATCH 5/5] separated some of the logic for jacobian calculation into separate function so it can be accessed from python. --- README.md | 0 example_scripts/dicom2nifti_example.sh | 0 example_scripts/run_SM_fit_example.py | 0 example_scripts/validate_formulas_example.py | 0 linear_planar/__init__.py | 0 linear_planar/b02b0_2mm.cnf | 0 linear_planar/calc_jacobian.py | 28 ++++++++++---------- linear_planar/debias.py | 0 linear_planar/eddy_nlgc_mp.py | 0 linear_planar/gs.py | 0 linear_planar/mean_std_4th.py | 0 linear_planar/smt_lin_pla.py | 0 linear_planar/split_shells_and_types.py | 0 preprocessing.sh | 0 14 files changed, 14 insertions(+), 14 deletions(-) mode change 100644 => 100755 README.md mode change 100644 => 100755 example_scripts/dicom2nifti_example.sh mode change 100644 => 100755 example_scripts/run_SM_fit_example.py mode change 100644 => 100755 example_scripts/validate_formulas_example.py mode change 100644 => 100755 linear_planar/__init__.py mode change 100644 => 100755 linear_planar/b02b0_2mm.cnf mode change 100644 => 100755 linear_planar/calc_jacobian.py mode change 100644 => 100755 linear_planar/debias.py mode change 100644 => 100755 linear_planar/eddy_nlgc_mp.py mode change 100644 => 100755 linear_planar/gs.py mode change 100644 => 100755 linear_planar/mean_std_4th.py mode change 100644 => 100755 linear_planar/smt_lin_pla.py mode change 100644 => 100755 linear_planar/split_shells_and_types.py mode change 100644 => 100755 preprocessing.sh diff --git a/README.md b/README.md old mode 100644 new mode 100755 diff --git a/example_scripts/dicom2nifti_example.sh b/example_scripts/dicom2nifti_example.sh old mode 100644 new mode 100755 diff --git a/example_scripts/run_SM_fit_example.py b/example_scripts/run_SM_fit_example.py old mode 100644 new mode 100755 diff --git a/example_scripts/validate_formulas_example.py b/example_scripts/validate_formulas_example.py old mode 100644 new mode 100755 diff --git a/linear_planar/__init__.py b/linear_planar/__init__.py old mode 100644 new mode 100755 diff --git a/linear_planar/b02b0_2mm.cnf b/linear_planar/b02b0_2mm.cnf old mode 100644 new mode 100755 diff --git a/linear_planar/calc_jacobian.py b/linear_planar/calc_jacobian.py old mode 100644 new mode 100755 index 8cb4504..23f0666 --- a/linear_planar/calc_jacobian.py +++ b/linear_planar/calc_jacobian.py @@ -48,7 +48,19 @@ def evaluate_jacobian(F, pixdim): return jacdet -def main(): +def calc_jacobian(in_file: str, out_file: str): + # Read data and extract necessary information + field = nib.load(in_file) + field_data = field.get_fdata() + hd_pixdim = field.header.get('pixdim') + pixdim = hd_pixdim[1:4] + + jacdet = evaluate_jacobian(field_data, pixdim) + + nib.save(nib.Nifti1Image(jacdet, field.affine), out_file) + + +if __name__ == '__main__': # Load parser to read data from command line input parser = buildArgsParser() args = parser.parse_args() @@ -60,16 +72,4 @@ def main(): DATA = os.path.realpath(args.data) OUT = args.out - # Read data and extract necessary information - field = nib.load(DATA) - field_data = field.get_data() - hd_pixdim = field.header.get('pixdim') - pixdim = hd_pixdim[1:4] - - jacdet = evaluate_jacobian(field_data, pixdim) - - nib.save(nib.Nifti1Image(jacdet, field.affine), OUT) - - -if __name__ == '__main__': - main() + calc_jacobian(DATA, OUT) diff --git a/linear_planar/debias.py b/linear_planar/debias.py old mode 100644 new mode 100755 diff --git a/linear_planar/eddy_nlgc_mp.py b/linear_planar/eddy_nlgc_mp.py old mode 100644 new mode 100755 diff --git a/linear_planar/gs.py b/linear_planar/gs.py old mode 100644 new mode 100755 diff --git a/linear_planar/mean_std_4th.py b/linear_planar/mean_std_4th.py old mode 100644 new mode 100755 diff --git a/linear_planar/smt_lin_pla.py b/linear_planar/smt_lin_pla.py old mode 100644 new mode 100755 diff --git a/linear_planar/split_shells_and_types.py b/linear_planar/split_shells_and_types.py old mode 100644 new mode 100755 diff --git a/preprocessing.sh b/preprocessing.sh old mode 100644 new mode 100755