Fix DKI radial kurtosis calcualtion and update WMTi

CHANGELOG.rst

@@ -4,7 +4,29 @@ Changelog
 All notable changes to this project will be documented in this file or
 page
 
-`v1.0_RC6`_
+`v1.0-RC7`_
+-----------
+
+Feb 11, 2021
+
+**Added**:
+
+* Missing Docker figures in RTD documentation
+
+* None
+
+**Changed**
+
+* Added error mitigation when FBI cost function fails to converge to
+  a minimum cost
+* Updated WMTI calculation to follow DKE outputs
+
+**Removed**
+
+* Unnecessary WMTI calculations
+
+
+`v1.0-RC6`_
 -----------
 
 Dec 22, 2020
@@ -29,7 +51,8 @@ Dec 22, 2020
 
 * See above
 
-`v1.0_RC5`_
+
+`v1.0-RC5`_
 -----------
 
 Oct 26, 2020
@@ -58,7 +81,7 @@ Oct 26, 2020
 
 * None
 
-`v1.0_RC4`_
+`v1.0-RC4`_
 -----------
 
 Sep 22, 2020
@@ -78,7 +101,7 @@ Sep 22, 2020
 
 * None
 
-`v1.0_RC3`_
+`v1.0-RC3`_
 -----------
 
 Sep 21, 2020
@@ -100,7 +123,7 @@ Sep 21, 2020
 
 * None
 
-`v1.0_RC2`_
+`v1.0-RC2`_
 -----------
 
 Aug 25, 2020
@@ -122,7 +145,7 @@ Aug 25, 2020
 
 * None
 
-`v1.0_RC1`_
+`v1.0-RC1`_
 -----------
 
 Aug 19, 2020
@@ -243,7 +266,7 @@ Feb 26, 2020
 
 * CSF masking; feature failed to work consistently
 
-`v0.11-dev`_
+`dev-0.11`_
 ------------
 
 Dec 2, 2019
@@ -262,7 +285,7 @@ Dec 2, 2019
 
 * None
 
-`v0.1-dev`_
+`0.1-dev`_
 -----------
 
 Oct 22, 2019
@@ -272,15 +295,16 @@ Initial port of MATLAB code to Python. 200,000,000,000 BCE
 
 .. Links
 
-.. _v1.0_RC6: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC6
-.. _v1.0_RC5: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC5
-.. _v1.0_RC4: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC4
-.. _v1.0_RC3: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC3
-.. _v1.0_RC2: https://github.com/m-ama/PyDesigner/releases/tag/v1.0_RC2
-.. _v1.0_RC1: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC1
+.. _v1.0-RC7: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC7
+.. _v1.0-RC6: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC6
+.. _v1.0-RC5: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC5
+.. _v1.0-RC4: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC4
+.. _v1.0-RC3: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC3
+.. _v1.0-RC2: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC2
+.. _v1.0-RC1: https://github.com/m-ama/PyDesigner/releases/tag/v1.0-RC1
 .. _v0.32: https://github.com/m-ama/PyDesigner/releases/tag/v0.32
 .. _v0.31: https://github.com/m-ama/PyDesigner/releases/tag/v0.31
 .. _v0.3: https://github.com/m-ama/PyDesigner/releases/tag/v0.3
 .. _v0.2 [The Cupid Release]: https://github.com/m-ama/PyDesigner/releases/tag/v0.2
-.. _v0.11-dev: https://github.com/m-ama/PyDesigner/releases/tag/dev-0.11
-.. _v0.2-dev: https://github.com/m-ama/PyDesigner/releases/tag/0.1-dev
+.. _dev-0.11: https://github.com/m-ama/PyDesigner/releases/tag/dev-0.11
+.. _0.1-dev: https://github.com/m-ama/PyDesigner/releases/tag/0.1-dev

designer/Untitled.ipynb

@@ -0,0 +1,91 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from designer.preprocessing import util, mrinfoutil"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bvec_new = '/media/sid/Secondary/Datasets/IAM_HiFI/out_test/working.mif'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 6,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "[0, 1000, 2000, 6000]\n"
+     ]
+    },
+    {
+     "data": {
+      "text/plain": [
+       "6000"
+      ]
+     },
+     "execution_count": 6,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "mrinfoutil.max_shell(bvec_new)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 7,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "True"
+      ]
+     },
+     "execution_count": 7,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "mrinfoutil.is_fullsphere(bvec_new)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.7.7"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}

designer/fitting/dwipy.py

@@ -674,6 +674,7 @@ def dkiTensorParams(self, v1, dt):
         akc = self.kurtosisCoeff(dt, dirs)
         ak = np.mean(akc)
         dirs = self.radialSampling(v1, dirSample)
+        akc = self.kurtosisCoeff(dt, dirs)
         rk = np.mean(akc)
         W_F = np.sqrt(dt[6]**2 + \
                       dt[16]**2 + \
@@ -1266,12 +1267,13 @@ def costCalculator(grid, BT, GT, b0, IMG, iDT, iaDT, zeta, shB, Pl0, g2l_fa_R_b,
                 raise Exception('Grid needs to be a flattened 1D vector')
             ndir = [len(x) for x in GT]
             cost_fn = np.zeros_like(grid)
-            for idx, awf in np.ndenumerate(grid):
-                for b in range(0, len(BT)):
-                    Se = (b0 * np.exp((-BT[b] * (1-awf)**-1) * np.diag((GT[b].dot((iDT - (awf**3 * zeta**-2) * iaDT).dot(GT[b].T)))))) * (1 - awf) # Eq. 3 FBWM paper
-                    Sa = (2*np.pi*b0*zeta*np.sqrt(np.pi/BT[b])) * (shB[b].dot((Pl0 * g2l_fa_R_b[b,idx,:][0]*clm))) # Eq. 4 FBM paper
-                    cost_fn[idx] = cost_fn[idx] + ndir[b]**-1 * np.sum((IMG[b] - Se.real - Sa.real)**2)
-                cost_fn[idx] = b0**-1 * np.sqrt(len(BT)**-1 * cost_fn[idx]) # Eq. 21 FBWM paper
+            with np.errstate(all='ignore'):
+                for idx, awf in np.ndenumerate(grid):
+                    for b in range(0, len(BT)):
+                        Se = (b0 * np.exp((-BT[b] * (1-awf)**-1) * np.diag((GT[b].dot((iDT - (awf**3 * zeta**-2) * iaDT).dot(GT[b].T)))))) * (1 - awf) # Eq. 3 FBWM paper
+                        Sa = (2*np.pi*b0*zeta*np.sqrt(np.pi/BT[b])) * (shB[b].dot((Pl0 * g2l_fa_R_b[b,idx,:][0]*clm))) # Eq. 4 FBM paper
+                        cost_fn[idx] = cost_fn[idx] + ndir[b]**-1 * np.sum((IMG[b] - Se.real - Sa.real)**2)
+                    cost_fn[idx] = b0**-1 * np.sqrt(len(BT)**-1 * cost_fn[idx]) # Eq. 21 FBWM paper
             return cost_fn
 
         def fbi_helper(dwi, b0, B, H, Pl0, gl, rectify=True,
@@ -1663,7 +1665,7 @@ def wmtihelper(dt, dir, adc, akc, awf, adc2dt):
             with np.errstate(invalid='ignore'):
                 akc[akc < minZero] = minZero 
             try:
-                # Eigenvalue decomposition of De
+                # Eigenvalue decomposition of De(extra-axonal)
                 De = np.multiply(
                     adc,
                     1 + np.sqrt(
@@ -1675,18 +1677,16 @@ def wmtihelper(dt, dir, adc, akc, awf, adc2dt):
                 eigval = np.sort(eigval)[::-1]
                 eas_ad = eigval[0]
                 eas_rd = 0.5 * (eigval[1] + eigval[2])
-                eas_md = np.add(eas_ad, (2 * eas_rd)) / 3
                 try:
                     eas_tort = eas_ad / eas_rd
                 except:
                     eas_tort = minZero
             except:
                 eas_ad = minZero
                 eas_rd = minZero
-                eas_md = minZero
                 eas_tort = minZero
             try:
-                # Eigenvalue decomposition of Da
+                # Eigenvalue decomposition of Da (intra-axonal)
                 Di = np.multiply(
                     adc,
                     1 - np.sqrt(
@@ -1696,20 +1696,11 @@ def wmtihelper(dt, dir, adc, akc, awf, adc2dt):
                 DTi = np.reshape(DTi, (3, 3), order='F')
                 eigval = sla.eigh(DTi, eigvals_only=True)
                 eigval = np.sort(eigval)[::-1]
-                ias_ad = eigval[0]
-                ias_rd = 0.5 * (eigval[1] + eigval[2])
-                ias_da = np.add(ias_ad, (2 * ias_rd))
+                ias_da = np.sum(eigval)
                 np.seterr(invalid='raise')
-                try:
-                    ias_tort = ias_ad / ias_rd
-                except:
-                    ias_tort = minZero
             except:
-                ias_ad = minZero
-                ias_rd = minZero
                 ias_da = minZero
-                ias_tort = minZero
-            return eas_ad, eas_rd, eas_md, eas_tort, ias_ad, ias_rd, ias_da, ias_tort
+            return eas_ad, eas_rd, eas_tort, ias_da
         dir = dwidirs.dirs10000
         nvox = self.dt.shape[1]
         N = dir.shape[0]
@@ -1751,7 +1742,7 @@ def wmtihelper(dt, dir, adc, akc, awf, adc2dt):
                       bar_format='{desc}: [{percentage:0.0f}%]',
                       unit='vox',
                       ncols=tqdmWidth)
-        eas_ad, eas_rd, eas_md, eas_tort, ias_ad, ias_rd, ias_da, ias_tort = zip(*Parallel(
+        eas_ad, eas_rd, eas_tort, ias_da = zip(*Parallel(
             n_jobs=self.workers, prefer='processes')(
             delayed(wmtihelper)(self.dt[:, i],
                                 dirs,
@@ -1768,7 +1759,7 @@ def wmtihelper(dt, dir, adc, akc, awf, adc2dt):
         ias_rd = vectorize(np.array(ias_rd), self.mask)
         ias_da = vectorize(np.array(ias_da), self.mask)
         ias_tort = vectorize(np.array(ias_tort), self.mask)
-        return awf, eas_ad, eas_rd, eas_md, eas_tort, ias_ad, ias_rd, ias_da, ias_tort
+        return awf, eas_ad, eas_rd, eas_tort, ias_da
 
     def findViols(self, c=[0, 1, 0]):
         """

designer/info.py

@@ -10,7 +10,7 @@
     )
 )
 __packagename__ = 'PyDesigner'
-__version__='v1.0-RC6'
+__version__='v1.0-RC7'
 __author__ = 'PyDesigner developers'
 __copyright__ = 'Copyright 2020, PyDesigner developers, MUSC Advanced Image Analysis (MAMA)'
 __credits__ = [