Lets assumed that $r$ is a vector that contains the residuals $r_i$ of the free water DTI model for all diffusion-weighted samples of an image voxel $s_i$. These residuals can be easly computed by subtracting each residual by the free water DTI model:
$$r_i = s_i - (1-x_8)\exp\left [ \sum_{j=1}^{7} W_{ij} x_j \right ] - x_8\exp\left [ \sum_{j=1}^{7} W_{ij} a_j \right ]$$where $x$ contain the six independent elements of the tissue diffusion tensor, the logarithm of the signal intensity when no diffusion-weighting is applied $S_0$, and the water volume fraction $f$:
$$x = \begin{bmatrix} D_{xx} & D_{xy} & D_{yy} & D_{xz} & D_{yz} & D_{zz} & \ln(S_0) & f \end{bmatrix} $$,
Based on the diffusion coefficient of free water at $37^oC$ ($D_{iso} = 3 \times 10 ^{-3} mm^{2}s^{-1} $), $a$ is defined as:
$$a = \begin{bmatrix} D_{iso} & 0 & D_{iso} & 0 & 0 & D_{iso} & \ln(S_0) \end{bmatrix} $$$W$ is computed from the diffusion $b_i$ values and diffusion gradient directions $g_i$ used to acquired the diffusion-weighted samples $s_i$:
$$ W = \begin{bmatrix} -b_1 g_{1x}^{2} & -2b_1 g_{1x}g_{1y} & -b_1 g_{2x}^{2} & -2b_1 g_{1x}g_{1z} & -2b_1 g_{1y}g_{1z} & -b_1 g_{1z}^{2} & 1\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ -b_m g_{mx}^{2} & -2b_m g_{mx}g_{my} & -b_m g_{mx}^{2} & -2b_m g_{mx}g_{mz} & -2b_m g_{my}g_{mz} & -b_m g_{mz}^{2} & 1 \end{bmatrix} $$.
The Jacobian matrix from scipy functions scipy.optimize.leastsq and scipy.optimize.least_squares are computed using:
The analytical solution of the partial derivatives are:
$$\left\{\begin{matrix} \frac{\partial r_i }{\partial x_j}= -(1-x_8)\exp\left [ \sum_{l=1}^{7} W_{il} x_l \right ]W_{ij}, &1\leq j\leq 6 \\ \frac{\partial r_i }{\partial x_j}= -(1-x_8)\exp\left [ \sum_{l=1}^{7} W_{il} x_l \right ]W_{ij}-x_8\exp\left [ \sum_{l=1}^{7} W_{il} a_l \right ]W_{ij}, & j=7 \\ \frac{\partial r_i }{\partial x_j}= \exp\left [ \sum_{l=1}^{7} W_{il} x_l \right ]-\exp\left [ \sum_{l=1}^{7} W_{il} a_l \right ], & j=8 \\ \\ \end{matrix}\right.$$When the free water volume fraction $f$ is converted to $f_t = \arcsin (2f-1) + \pi / 2$ to insure its convergence within a plausible range, $f$ in the last element of vector x has to be replaced by $f_t$ and the residual function converted to:
$$r_i = s_i - [1- 0.5 (1 + \sin(x_8 - \pi/2))]\exp\left [ \sum_{j=1}^{7} W_{ij} x_j \right ] - 0.5 (1 + \sin(x_8 - \pi/2))\exp\left [ \sum_{j=1}^{7} W_{ij} a_j \right ]$$In this case the Jacobian matrix solution is:
$$\left\{\begin{matrix} \frac{\partial r_i }{\partial x_j}= -(1-x_8)\exp\left [ \sum_{l=1}^{7} W_{il} x_l \right ]W_{ij}, &1\leq j\leq 6 \\ \frac{\partial r_i }{\partial x_j}= -(1-x_8)\exp\left [ \sum_{l=1}^{7} W_{il} x_l \right ]W_{ij}-x_8\exp\left [ \sum_{l=1}^{7} W_{il} a_l \right ]W_{ij}, & j=7 \\ \frac{\partial r_i }{\partial x_j}= 0.5*\cos(x_8-\pi/2)\exp\left [ \sum_{l=1}^{7} W_{il} x_l \right ]-0.5*\cos(x_8-\pi/2)\exp\left [ \sum_{l=1}^{7} W_{il} a_l \right ], & j=8 \\ \\ \end{matrix}\right.$$
In [ ]: